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Introduction 


There  is  currently  no  cure  or  dependable  prevention  method  for  breast  cancer,  and  early 
detection  by  breast  imaging  provides  the  current  best  practice  for  providing  effective  treatment  and 
reducing  mortality.  While  conventional  x-ray  mammography  reduces  breast  cancer  mortality  by 
about  30%  in  women  over  50  and  by  17%  in  women  ages  40-50,  x-ray  methods  have  limited 
sensitivity  and  specificity  for  breast  cancer  detection  and  diagnosis.  Many  women  have  radiodense 
breast  tissue,  and  in  monitoring  breast  cancer  survivors  for  recurrence  scar  tissue  can  make  x-ray 
image  analysis  very  difficult.  Overall,  roughly  75%  of  lesions  detected  by  conventional 
mammography  are  benign,  leading  to  a  large  number  of  unnecessary  biopsies;  in  addition,  x-ray 
techniques  give  15-20%  false  negative  results  when  cancerous  lesions  are  in  fact  present  [1]. 
These  limitations  in  conventional  mammography  have  given  rise  to  an  urgent  need  to  develop  novel 
imaging  technologies  for  more  effective  early  breast  cancer  detection. 

Among  potential  alternative  breast  imaging  technologies,  Positron  Emission  Tomography 
(PET)  offers  a  unique  combination  of  demonstrated  high  sensitivity  and  high  specificity  for  breast 
cancer  detection,  plus  the  potential  for  providing  new  in  vivo  imaging  capabilities  building  on  our 
improving  understanding  of  the  molecular  biology  of  breast  cancer.[2]  PET  is  a  functional 
imaging  modality,  and  as  such  is  complimentary  to  the  anatomical  imaging  provided  by  x-ray, 
ultrasound,  and  MRI  techniques.  In  PET,  the  patient  is  injected  with  a  short-lived  positron- 
emitting  radiotracer  which  is  preferentially  taken  up  by  malignant  tissues;  since  positron-emitting 
isotopes  of  Carbon,  Oxygen,  Nitrogen  and  Fluorine  (which  serves  as  a  chemical  analog  of 
Hydrogen)  exist,  positron-emitting  forms  of  virtually  any  biomolecule  may  be  produced,  at  least  in 
principle.  While  the  great  majority  of  clinical  PET  studies  of  breast  cancer  have  been  performed 
with  a  radioactive  form  of  glucose  (FDG  =  18F-Fluoro-Deoxy-Glucose)  [3-8],  a  number  of 
pioneering  studies  have  used  radioactive  taxol,  estrogens,  and  various  other  ligands  [9-11]. 

PET  breast  imaging  can  exploit  a  panoply  of  detection  instruments  as  well  as  a  variety  of 
radiotracers;  in  addition,  nuclear  medicine  provides  a  number  of  additional  lower-energy  planar 
imaging  and  SPECT  (Single  Photon  Emission  Computed  Tomography)  modalities  and  radiotracers 
under  the  general  appelation  “scintimammography”.[12]  In  addition  to  its  great  potential  for 
radiotracer  diversity,  PET  is  distinguished  from  scintimammography  by  higher  spatial  resolution 
and  lower  scatter  backgrounds,  which  should  lead  to  higher  sensitivity  and  specificity  for  small, 
early  lesions.  Within  the  last  few  years,  several  research  groups  have  begun  to  explore  the 
potential  of  dedicated  PET  breast  imaging  systems,  whose  greater  proximity  to  the  breast  should 
lead  to  improved  counting  statistics  and  higher-quality  images.  The  reduced  detector  size  for  a 
dedicated  system  typically  allows  for  instrumentation  choices  which  would  be  quite  expensive  in 
large  imaging  systems,  again  giving  the  potential  for  higher-resolution  and  higher-quality  images. 
At  the  time  of  this  writing,  we  are  aware  of  5  groups  (plus  our  own)  who  have  proposed  or  are 
constructing  dedicated  PET  breast  imaging  systems,  but  no  clinical  results  have  thus  far  been 
presented  or  published  by  any  group  [13-17].  The  objective  of  our  ongoing  research  effort  is  to 
construct  and  characterize  components  of  a  dedicated  PET  breast  imaging  system  which  we  believe 
will  have  technical  advantages  over  its  rivals,  where  these  advantages  should  lead  to  improved 
image  quality  and  quantitative  accuracy  at  lesser  cost,  as  well  as  greater  sensitivity  and  greater 
accuracy  near  the  physical  edge  of  the  detector  (for  imaging  axillae  and  near  the  chest  wall). 

Our  two-year  research  effort  is  focussed  on  the  design  and  fabrication  of  high-resolution 
PET  mammography  modules  according  to  our  novel  design,  and  on  the  evaluation  and 
optimization  of  their  anticipated  performance  as  part  of  a  future  dedicated  PET  mammograph. 
Within  the  past  year  we  have  constructed  and  tested  a  succession  of  prototype  modules,  and  have 
evaluated  their  performance  in  response  to  sources  with  known  characteristics  and  geometry.  The 
most  important  features  of  our  module  design  and  the  principal  results  of  key  performance 
assessments  will  be  detailed  in  the  next  section,  followed  by  a  discussion  of  our  plans  for  the 
further  assessment  of  these  and  successor  modules  within  the  coming  year. 
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Experimental  Methods,  Assumptions  and  Procedures: 

Our  detector  design  is  distinguished  by  its  geometry  and  by  its  readout  method,  which 
combine  to  give  improved  high-resolution  3-D  imaging  capability  at  reduced  cost  in  comparison  to 
rival  techniques.  Conventional  PET  detector  module  designs  use  scintillator  crystals  to  serve  two 
purposes:  to  convert  the  energy  of  incident  gamma  rays  to  visible  light,  and  to  guide  this  visible 
light  to  one  or  more  photosensors  which  are  optically  coupled  to  the  crystals.  In  a  conventional 
detector,  one  uses  an  array  of  small  scintillator  crystals  and  various  optics/electronics  designed  to 
ascertain  which  crystal  converted  the  gamma  ray,  with  little  or  no  information  on  where  within  the 
crystal  such  conversion  occurred.  Typical  scintillator  crystal  geometries  are  perhaps  3mm  x  3mm 
x  20mm,  where  the  last  dimension  is  dictated  by  the  need  for  good  conversion  efficiency  for  the 
relatively  high-energy  511keV  gammas.  A  pair  of  such  scintillator  arrays  is  needed  for  PET 
imaging,  since  the  radiotracer  source  distribution  is  tomographically  reconstructed  from  lines-of- 
response  connecting  pairs  of  coincident  gamma  ray  interactions  (where  these  gamma  rays  have 
arisen  from  positron  annihilation  and  the  emission  of  back-to-back  gammas). 

Since  the  range  of  a  secondaiy  electron  produced  in  the  photocapture  reaction  of  a  gamma 
ray  within  a  dense  scintillator  crystal  is  much  less  than  a  millimeter,  considerable  spatial 
information  is  lost  if  the  position  of  interaction  within  a  given  scintillator  crystal  is  not  measured 
(as  it  is  not,  in  conventional  detectors).  Most  importantly,  gamma  rays  which  are  incident  on  the 
scintillator  array  at  oblique  angles  have  significant  errors  in  their  coincident  lines-of-response,  an 
effect  which  is  usually  termed  parallax  error  or  depth-of-interaction  ambiguity.  This  effect 
becomes  particularly  pronounced  in  the  close-proximity  geometry  of  a  dedicated  PET  breast 
imaging  system,  most  especially  if  one  seeks  to  use  coincident  events  with  fairly  large  oblique 
angles  of  incidence  in  order  to  obtain  significant  spatial  resolution  along  the  axis  connecting  two 
opposed  detector  arrays.  Our  competitors’  dedicated  PET  breast  imaging  designs,  therefore,  either 
reduce  the  depth  of  the  scintillator  crystals  to  much  less  than  20mm  and  thus  have  lesser  efficiency, 
or  else  reject  events  at  large  incident  angles  in  what  amounts  to  planar  rather  than  tomographic 
imaging.  The  former  strategy  reduces  event  statistics  and  thereby  reduces  spatisl  resolution  given 
reasonable  limits  on  the  amount  of  radiotracer  injected  in  the  patient,  while  the  latter  strategy 
reduces  image  signal-to-noise  by  superimposing  overlying  tissue  on  the  region  of  interest. 

Our  design  avoids  depth-of-interaction  ambiguity  by  measuring  gamma-ray  interaction 
positions  before  scintillation  light  has  spread  far;  this  is  accomplished  by  configuring  our 
scintillator  as  a  stack  of  thin,  flat  sheets  (2.5mm  x  112mm  x  1 12mm)  interspersed  with  ribbons  of 
wavelength-shifting  (fluorescent)  optical  fibers.  These  fibers  serve  as  photonic  sensors,  taking 
scintillation  light  as  input  and  providing  optically-piped  fluorescent  light  as  output.  This 
wavelength-shifted  tight  is  then  guided  to  position-sensitive  photosensors  for  readout,  thereby 
measuring  the  position  of  the  gamma-ray  interaction  and  scintillation  tight  production.  Unlike 
transparent  fiber  optics,  fluorescent  fiber  optics  can  accept  input  light  through  their  sides  and  pipe  it 
to  their  ends,  so  that  a  “ribbon”  of  parallel  fibers  can  serve  as  a  incident  tight  position  sensor  in  1 
dimension.  We  use  perpendicular  ribbons  of  wavelength-shifting  fibers  on  opposite  sides  of  our 
thin  scintillator  sheets  to  determine  gamma-ray  interaction  positions  in  two  dimensions,  with  the 
third  interaction  coordinate  determined  up  to  the  crystal  thickness  within  a  stack.  For  a  polished 
scintillator  crystal,  only  scintillation  tight  travelling  nearly  normal  to  flat  surfaces  may  emerge  from 
the  crystal  without  undergoing  total  internal  reflection,  due  to  the  high  refractive  index  of  the 
crystal  relative  to  its  surroundings.  The  optical  principles  underlying  our  design  are  further 
detailed  in  attached  Appendix  I:  “Development  of  a  High-Resolution  PET  Detector  using  LSO  and 
Wavelength-Shifting  Fibers”,  from  the  Proceedings  of  the  1995  IEEE  Conference  on  Medical 
Imaging.  The  work  described  in  this  Appendix  directed  preceded  our  USAMARC-funded  effort, 
and  culminated  in  our  demonstration  of  very  high  spatial  resolution  for  small  LSO  crystals  coupled 
to  wavelength-shifting  fibers. 
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The  choice  of  optimal  materials  and  methods  for  the  construction  and  readout  of  high- 
resolution  PET  mammography  modules  with  wavelength-shifting  fiber  readout  has  been  a  major 
part  of  this  year’s  effort,  and  the  basis  for  our  design  choices  will  be  discussed  in  the  “Results” 
section  which  follows.  We  have  extensively  tested  the  available  design  options  at  each  step  in  the 
optical  and  readout  chain.  This  testing  has  been  carried  out  with  a  nearly  point-like  positron  source 
(Na-22)  which  was  collimated  for  testing  module  components  and  single  modules,  and  then  used 
for  coincidence  testing  of  pairs  of  prototype  modules.  The  critical  performance  parameters 
measured  included  system  spatial  resolution,  energy  resolution,  timing  resolution,  efficiency,  light 
yield  (which  correlates  with  the  robustness  and  efficiency  of  the  fiber  coordinate  measurement), 
and  rate  capability.  We  have  used  a  50  microcurie  source  for  most  tests,  where  this  is 
approximately  the  total  activity  expected  within  a  breast  during  a  PET  mammograph  scan. 

One  of  the  most  important  design  choices  comes  at  the  beginning  of  the  optical  chain,  with 
the  choice  of  scintillator.  In  our  proposal  we  planned  to  use  LSO  or  Lutetium  Oxyorthosilicate, 
and  we  have  been  able  to  obtain  lOOcc’s  of  this  material  for  the  construction  of  test  modules. 
Unfortunately,  this  proprietary  material  has  until  very  recently  (just  this  month)  been  unavailable 
outside  of  CTI/Siemen’s  internal  development  program,  so  that  we  have  switched  to  a  less 
expensive  and  more  readily  available  scintillator  for  our  initial  work.  We  selected  CsI(Na), 
Sodium-doped  Cesium  Iodide,  based  upon  its  high  light  yield,  ease  of  working,  and  relatively  low 
cost;  as  will  be  discussed  below  this  choice  has  placed  some  limits  on  our  achievable  spatial 
resolution  and  timing/rate  capability,  but  it  has  allowed  us  to  construct  many  prototypes  rapidly 
with  which  to  explore  design  issues  which  are  common  to  CsI(Na)  or  LSO-based  devices.  Given 
the  order-of-magnitude  lower  cost  of  CsI(Na)  relative  to  LSO,  there  is  significant  motivation  for 
obtaining  good  performance  with  a  CsI(Na)  device;  we  have  succeeded  in  obtaining  funding  for 
limited  clinical  application  of  a  larger-scale  successor  CsI(Na)  imaging  system  from  the 
Commonwealth  of  Massachusetts  Breast  Cancer  Research  Program  under  the  title  “Cost-Effective 
High-Resolution  Positron  Emission  Mammography.”  We  nonetheless  intend  to  continue 
development  of  modules  using  LSO  scintillator  as  well  as  CsI(Na)-based  modules,  and  we  have 
lOOcc’s  of  LSO  currently  in  hand  which  we  are  preparing  for  incorporation  into  two  single-layer 
modules  each  1 12mm  x  1 12mm  x  2mm  thickness.  Tests  on  these  LSO  modules  will  be  carried  out 
in  the  coming  year,  and  we  will  acquire  sufficient  LSO  for  two  3-layer  LSO  modules,  as  planned. 

Our  most  important  design  assumption  to  date  has  been  that  the  radiation  environment 
(source  strength  and  distribution,  background  rates,  etc.)  in  which  we  are  testing  our  prototypes  is 
comparable  to  that  which  will  be  experienced  in  an  actual  mammography  module  during  a  clinical 
PET  imaging  procedure.  Our  source  activity  estimates  are  based  upon  measurements  which  were 
performed  with  conventional  whole-body  PET  scanners  and  reported  in  the  literature,  combined 
with  estimates  of  activity  outside  the  field-of-view  (such  as  the  radiotracer  uptake  in  the  heart, 
particularly  for  FDG)  and  its  expected  contribution  to  backgrounds  by  scattering  into  our  detectors. 
These  are  rough  estimates  and  we  have  yet  to  carry  out  careful  analyses  of  optimal  screening  of 
activity  outside  the  field-of-view,  since  we  expect  the  most  useful  input  information  for  such  an 
analysis  to  come  from  experience  with  our  actual  modules  in  future  clinical  application.  Such 
future  measurements,  and  any  forms  of  human  use,  are  beyond  the  scope  of  our  research  program 
as  funded  by  the  USAMARC  and  detailed  in  our  statement  of  work,  but  planning  for  them  has  had 
a  significant  impact  on  our  detector  design.  For  example,  we  have  constructed  our  prototypes  to 
have  an  “active  edge”  and  an  “active  comer”  based  on  discussions  with  clinicians  who  emphasized 
the  importance  of  axillary  nodal  imaging  and  of  imaging  tissues  near  the  chest  wall.  We  have  made 
every  effort  to  make  our  device  geometry  compatible  with  breast  imaging  practice  as  carried  out  in 
x-ray  mammography  (although  we  anticipate  much  less  compression  to  accompany  our 
considerably  longer  scans),  and  we  have  gone  so  far  as  to  obtain  a  donated  mammography  stand 
which  we  are  outfitting  with  mechanical  fixtures  to  mount  our  prototype  detectors  and  their  clinical 
successors.  Our  intention  is  to  move  our  research  from  the  laboratory  to  the  clinic  at  the  earliest 
reasonable  date,  and  our  eyes  are  steady  on  this  future  transition. 
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Detector  Design  and  Construction: 

As  mentioned  earlier,  the  details  of  our  prototype  detector  design  have  evolved  through 
several  iterations  of  a  construction/test  cycle,  and  have  been  informed  by  a  variety  of  component- 
level  tests  and  computer  simulations.  The  experimental  findings  upon  which  we  have  based  our 
most  significant  design  decisions  are  detailed  below: 

1)  Crystal  geometry  and  surface  preparation: 

Our  crystal  geometry  of  112mm  x  112mm  x  2.5mm  was  based  upon  readout  through  8 
ribbons  each  on  the  x-  and  y-  surfaces  of  each  crystal,  where  each  ribbon  is  14mm  across.  This 
14mm  number  arose  from  measurements  on  early  production  units  of  the  Hamamatsu  R5900-L16 
position-sensitive  photomultipliers  which  we  use  to  read  out  the  fibers;  more  recent  units  have 
shown  good  performance  for  16mm  width  ribbons,  so  that  future  crystal  dimensions  of  128mm  x 
128mm  x  2.5mm  are  appropriate.  The  2.5mm  thickness  was  based  upon  a  target  resolution  of 
2.0mm  FWHM,  which  we  expect  to  be  a  good  match  to  the  statistical  limits  in  image  quality 
available  for  a  high-efficiency  PET  mammography  unit  for  anticipated  clinical  doses  and  scan 
times.  The  actual  resolution  obtained  with  CsI(Na)  is  2.5mm,  due  to  Compton  scatter  interactions 
within  the  large,  flat  scintillator  crystals  as  discussed  further  below.  Rather  than  decrease  the 
crystal  thickness  further,  which  would  increase  the  number  of  layers  required  for  an  efficient 
detector,  we  expect  to  keep  the  same  thickness  until  we  have  more  experience  with  clinical 
experience  with  successor  prototypes.  For  LSO-based  modules,  we  would  anticipate  a  crystal 
thickness  of  2.0mm  in  order  to  take  advantage  of  the  higher  density,  faster  timing,  and  higher 
photofraction  of  this  material  relative  to  CsI(Na);  in  both  cases,  a  full-thickness  module  is  expected 
to  incorporate  7  layers  of  crystal  read  out  through  8  fiber  ribbons. 

The  quality  of  the  scintillator  crystal  surface  finish  is  important  in  our  application  for  two 
reasons:  first,  because  we  require  a  clear  division  between  light  which  is  totally  internally  reflected 
within  the  crystal  and  that  light  which  emerges  from  the  crystal  into  the  fiber  ribbons,  and  second 
because  we  need  the  overall  crystal/fiber  laminated  stack  to  be  transparent  to  wavelength-shifted 
light.  We  require  such  transparency  because  in  our  device  we  have  separated  the  functions  of 
interaction  position  measurement  (through  the  fiber  readout)  from  gamma  ray  energy  measurement 
and  triggering/timing  (through  an  Anger  array  of  photomultipliers  which  view  each  stack). 
Because  both  the  scintillator  crystals  and  the  fibers  are  transparent  to  the  wavelength-shifted  light, 
photons  traversing  the  wide  range  of  optical  indices  within  the  stack  (fiber  core,  fiber  claddings, 
optical  coupling  agent,  scintillator  coating,  and  scintillator  crystal)  are  scattered  rather  than 
absorbed;  small  scratches  on  the  crystal  surface  from  imperfect  polishing,  however,  absorb  light 
and  can  cause  excessive  attenuation  of  wavelength-shifted  light  within  the  laminate.  In  our 
devices,  we  surround  the  laminated  stack  by  reflective  material  on  all  sides  except  that  side  with  the 
Anger  array,  and  we  are  able  to  obtain  good  energy  and  timing  resolution  as  well  as  good 
uniformity  of  response.  Because  CsI(Na)  does  not  readily  polish  to  an  optical  finish  and  because  it 
is  somewhat  hygroscopic,  considerable  effort  was  needed  in  order  to  obtain  an  adequate  surface 
finish  (using  silicon  oil  and  progressively  smaller  sizes  of  Aluminum  Oxide  powder  abrasives). 
Our  final  procedure  yielded  less  than  10-20%  variation  in  light  collection  from  one  side  of  the  stack 
versus  the  other,  which  is  readily  corrected  since  the  depth  of  interaction  in  the  stack  is  measured. 

The  effect  of  imperfect  surface  preparation  for  CsI(Na)  was  more  difficult  to  analyze,  and 
determination  of  the  role  of  scintillation  light  scattering  from  crystal  surface  imperfections  (rather 
than  totally  internally  reflecting)  required  a  great  deal  of  experimentation.  These  experiments  were 
complicated  by  the  fact  that  the  radioactive  source  we  were  using  to  assess  our  system’s  spatial 
resolution  turned  out  to  have  considerable  extension;  although  specified  to  have  active  dimensions 
of  1mm  x  1mm  x  1mm  at  the  time  of  its  purchase,  it  turned  out  to  be  a  disk  roughly  5mm  in 
diameter  and  perhaps  1mm  in  thickness.  Once  we  realized  that  we  had  been  imaging  our  sources’s 
extension  rather  than  measuring  our  device’s  resolution,  we  had  built  special  test  modules  with 
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lmm-deep  reflector-filled  grooves  at  2mm  intervals  across  each  crystal  surface  —  in  order  to  limit 
light  spreading  away  from  the  interaction  region  within  each  scintillator  crystal.  When  imaging  the 
smallest  dimension  of  our  disk-shaped  source,  we  then  found  no  difference  in  spatial  resolution 
between  the  grooved  modules  and  earlier  ungrooved  modules,  eliminating  the  hypothesis  that  our 
observed  2.5mm  limiting  resolution  was  due  to  light  spreading.  This  limiting  resolution  was  also 
insensitive  to  the  amount  of  light  collected  by  the  fibers,  which  was  a  tip-off  that  it  was  either  due 
to  source  extension  or  to  non-pointlike  light  production  in  the  scintillator  crystal. 

Finally,  we  were  able  to  determine  that  our  imaging  resolution  for  2.5mm  thick  plates  of 
CsI(Na)  is  limited  by  Compton  scattering  and  secondary  photon  capture  within  a  single  plate.  The 
true  photofraction  (ratio  of  photocapture  events  to  all  gamma-ray  interactions)  for  511  keV  gamma 
rays  in  CsI(Na)  is  about  15%,  while  the  apparent  photofraction  (fraction  of  events  within  the 
photopeak  for  51 1  keV  gammas  normally  incident  on  a  2.5mm  thick  plate)  is  more  than  25%.  The 
difference  arises  from  gammas  which  Compton  scatter  at  90  degrees,  followed  by  photocapture  of 
the  secondary  gamma  rays.  These  secondary  gammas  have  energies  near  250  keY  and  a  range  of 
about  7  mm  in  CsI(Na);  computer  simulations  show  a  limiting  resolution  due  to  this  effect  close  to 
our  observed  2.5mm  FWHM,  for  two  detectors  in  coincidence.  For  a  CsI(Na)  device,  this  effect 
can  be  somewhat  lessened  with  thinner  crystals,  which  increases  the  fraction  of  “true” 
photocaptures  within  “apparent”  photocaptures  (by  decreasing  the  angular  acceptance  for 
secondaries  to  convert  in  the  same  crystal)  but  this  either  increases  the  number  of  detector  layers 
required  or  decreases  the  overall  efficiency.  For  LSO  devices  the  effect  is  much  less  pronounced 
because  of  the  higher  true  photofraction  for  LSO  (30%)  and  the  factor  of  2  shorter  attenuation 
length  for  the  secondary  gammas  in  this  denser  material  relative  to  CsI(Na).  In  fact,  it  was  because 
our  earlier  tests  on  fiber  readout  of  LSO  had  shown  such  good  spatial  resolution  that  we  were 
somewhat  surprised  by  the  magnitude  of  this  effect  in  CsI(Na). 

2)  Optical  coupling  of  crystals  to  fibers: 

Our  second  round  of  prototype  modules  (the  first  had  attenuation  difficulties  due  to 
imperfect  crystal  surface  treatment)  yielded  a  significantly  lower  fiber  light  yield  relative  to  the  first, 
for  reasons  which  were  not  understood  at  first.  After  some  experimentation,  we  discovered  that 
the  effect  was  due  to  a  relatively  high-index  optical  epoxy  which  we  used  to  coat  the  hygroscopic 
crystals;  once  we  determined  that  the  crystals  did  not  degrade  with  ambient  humidity  over  a  time 
scale  of  1-2  weeks  (the  time  between  cutting/polishing  and  sealing  them  into  laminated  modules) 
we  omitted  the  coating  step.  Although  the  crystals  did  not  degrade,  the  greater  optical  index 
mismatch  without  the  coating  trapped  more  light  within  the  crystal  and  decreased  our  fiber  light 
yield.  In  the  third  round  of  prototypes,  the  light  yield  was  restored  and  slightly  improved  by 
introducing  a  crystal  coating  layer  with  refractive  index  near  the  geometric  mean  between  that  of  the 
crystal  and  that  of  the  fiber  outer  cladding;  this  thin  layer  then  serves  as  an  antireflection  coating. 
There  is  still  sufficient  light  trapped  within  each  crystal  to  be  collected  at  the  perimeter  of  each 
crystal  within  a  stack,  however,  as  will  be  described  in  the  module  assembly  section  below. 

3)  Fiber  selection  and  preparation: 

We  experimented  with  both  round  and  square  cross-section  fibers,  with  different  fluor 
formulations  and  concentrations,  and  with  diameters  ranging  from  0.5mm  to  1.0mm.  Round 
fibers  gave  nearly  a  factor  of  two  better  light  yield  relative  to  square  fibers,  which  we  traced  to 
much  smaller  attenuation  losses  along  the  fibers  with  rouind  cross-sections.  This  arises  because 
the  square  fibers  do  not  have  a  precisely  square  cross-section,  leading  to  optical  ray  paths  which 
encounter  the  core/cladding  interface  at  more  than  the  critical  glancing  incidence  angle.  This  effect 
may  be  evidenced  by  shining  a  flashlight  on  the  side  of  a  fiber  ribbon  in  a  darkened  room  -  this 
flashlight  test  also  proved  useful  when  examining  detector  modules  for  any  fiber  damaage  such  as 
nicked  cladding.  Similarly,  1mm  fibers  had  significantly  higher  light  yield  than  0.5mm  fibers,  in 
this  case  due  to  self-absorption  to  the  limited  Stokes  shift  between  absorption  and  emission. 
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combined  with  the  higher  fluor  concentrations  required  for  efficienct  fluorescent  in  0.5mm  fibers. 
Finally,  BBQ-doped  fibers  did  not  yield  a  significantly  higher  light  yield  than  our  standard  BCF- 
91A  formulation  when  coupled  to  CsI(Na)  crystals.  We  have  used  multi-clad  fibers  in  all 
prototypes  (which  earlier  tests  showed  superior  to  single-clad  fibers)  and  succeeded  in  obtaining 
50%  light  yield  improvements  by  optically  coupling  aluminized  mylar  mirrors  to  the  ends  of 
polished  fiber  ribbons.  This  last  was  used  to  obtain  two  active  edges  and  an  active  corner  away 
from  the  readout  ends  of  the  fibers. 

41  Multilayer  stack  assembly  and  fixing: 

Multilayer  crystal/fiber  stacks  were  built  up  from  acrylic  “frames”  3.5mm  thick,  with 
cutaway  regions  on  one  side  to  guide  a  1mm  thickness  fiber  ribbon,  and  a  square  opening  in  the 
center  to  position  the  112mm  x  112mm  coated  crystal.  Fiber  ribbons  were  first  glued  to  frames, 
then  a  stack  was  built  up  by  placing  coated  crystals  onto  fibers  within  each  frame,  with  the  axes  of 
the  fiber  ribbons  alternating  between  two  orientations  as  the  stack  was  built  up.  We  then  fixed  the 
assembled  stacks  by  casting  them  in  an  epoxy  resin  (Emerson  &  Cummings  potting  compound 
Stycast  1367).  We  devoted  considerable  effort  to  developing  procedures  which  allowed  the 
potting  compound  to  fill  all  interstices  without  generating  air  gaps  or  bubbles,  and  we  succeeded  in 
casting  nearly  defect-free  7-layer  stacks.  For  each  stack,  the  top  layer  contained  fibers  only,  so 
that  a  stack  with  7  crystal  layers  has  8  fiber  layers;  this  is  required  in  order  to  obtain  x-  and  y- 
coordinate  information  from  all  crystals  in  the  stack.  Wavelength-shifting  fibers  were  also  coupled 
to  the  perimeter  of  each  crystal  in  a  stack,  with  these  fibers  read  out  separately  to  determine  the 
depth-of-interaction  within  a  stack  (these  fibers  detected  the  tight  which  was  trapped  in  a  crystal 
layer  and  guided  to  the  perimeter  through  total  internal  reflection).  We  obtained  a  tight  yield  of 
about  10  photoelectrons  for  0.5mm  diameter  perimeter  fibers,  which  is  more  than  sufficient  to 
determine  depth-of-interaction  to  within  the  crystal  thickness  of  2.5mm.  The  top  of  each  module 
was  epoxied  to  a  5cm-thick  “tight  mixer”  which  was  later  coupled  to  four  2”  diameter 
photomultipliers  in  an  Anger  array.  The  bottom  of  each  module  was  optically  coupled  to  reflecting 
material  after  assembly,  as  were  the  edges  of  the  tight  mixer  and  those  regions  of  the  tight  mixer 
not  covered  by  Anger  PMTs.  Fiber  ends  which  were  not  read  out  were  polished  and  optically 
coupled  to  aluminized  mylar  mirrors. 

5)  Fiber  ribbon  routing,  multiplexing,  and  readout  coupling: 

To  date  we  have  been  reading  out  only  single-layer  modules  (i.e.  a  single  crystal  layer 
coupled  to  8  x  14mm  width  x-ribbons  and  8  x  14mm  width  y-ribbons).  Initial  testing  was 
performed  with  no  readout  multiplexing,  reading  out  just  4  x-ribbons  and  4  y-ribbons  through  8 
multianode  photomultipliers.  This  was  followed  by  multiplexed  readout  of  two  ribbons  per 
multianode  photomultiplier,  where  a  rough  coordinate  from  the  Anger  PMT  array  was  used  to 
demultiplex  between  the  two  ribbons  on  each  fiber  readout  PMT.  For  our  initial  Anger  optics  it 
was  two  ribbons  4  x  14mm  apart  were  very  easily  distinguished,  as  were  two  ribbons  2  x  14mm 
apart  near  the  center  of  the  module.  However,  the  edgemost  ribbon  and  that  2  x  14mm  further 
toward  the  module  center  had  a  small  amount  of  overlap;  further  improvements  in  the  Anger  optics 
would  thus  be  required  for  four-fold  lateral  multiplexing  (i.e.  reading  out  8  ribbons  with  2 
position-sensitive  photomultipliers).  There  is  sufficient  photocathode  area  on  our  multianode 
photomultipliers  to  permit  both  4-fold  lateral  multiplexing  and  4-fold  multiplexing  in  depth  when 
using  1mm  diameter  fibers,  but  now  it  appears  such  heavy  multiplexing  will  not  be  needed  because 
of  the  new  multianode  photomultiplier  options  discussed  in  the  next  section.  Given  the  $1000 
price  of  each  multianode  photomultiplier  it  is  important  to  keep  their  numbers  down,  but  heavy 
multiplexing  decreases  system  rate  capability  by  introducing  pile-up  and  demultiplexing  errors  at 
high  event  rates.  This  is  much  more  the  case  for  CsI(Na),  with  its  700ns  decay  time,  than  for  LSO 
with  its  70ns  decay  time,  but  with  the  new  multianode  photomultipliers  we  can  greatly  reduce  the 
optical  multiplexing  while  keeping  PMT  counts  down. 
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In  all  cases,  fibers  were  optically  coupled  to  multianode  PMTs  by  inserting  fiber  ribbons 
into  custom-built  vises,  polishing  the  fiber  ends,  then  optically  and  mechanically  coupling  to  the 
multianode  PMTs  by  attaching  the  vises  to  custom  PMT  mounting  mechanics.  We  found  that 
coupling  fibers  to  PMT  photocathodes  with  silicon  optical  coupling  compound  was  sufficient  to 
reproducibly  produce  a  uniform  and  high-efficiency  coupling;  this  method  was  conducive  to  the 
testing  of  alternative  fiber  multiplexing  schemes  by  easily  making  and  breaking  optical  contacts. 

6)  Mnltianode  photomultipliers  and  associated  electronics: 

We  used  Hamamatsu  R5900-L16  multianode  photomultipliers  to  read  out  the  fibers,  where 
each  multianode  PMT  contained  16  separate  photocathode  regions  each  16mm  x  1mm.  The  16 
anodes  were  resistively  coupled  in  a  chain  and  read  out  through  charge  division;  details  of  this 
readout  and  measurements  of  the  PMT  response  are  given  in  Appendix  II:  “Characterization  of  a 
new  Multianode  PMT  for  Low  Light-level  Optical  Fiber  Readout”,  from  the  Proceedings  of  the 
1996  IEEE  Nuclear  Science  Symposium  and  Medical  Imaging  Conference.  We  were  able  to  obtain 
spatial  resolution  of  less  than  1mm  FWHM  for  8  photoelectrons  incident  through  a  single  1mm 
fiber;  our  spatial  resolution  is  therefore  not  limited  by  our  PMT  readout,  but  rather  by  upstream 
optics.  The  two  charge-division  outputs  (whose  ratio  give  the  transverse  coordinate  for  light 
incident  through  a  fiber  ribbon,  and  whose  sum  gives  the  total  fiber  light  yield)  were  read  out 
through  custom-built  4-PMT  bases  and  coupled  to  xlO  charge-sensitive  preamplifiers  which  we 
also  constructed.  The  single  photoelectron  peak  from  the  multianode  photomultipliers  was  very 
pronounced,  simplifying  detector  gain  calibration  considerably.  Spatial  calibration  of  the  detector 
was  also  straightforward  when  imaging  a  point  source,  and  3-4mm  FWHM  resolution  was 
obtainable  with  no  calibration  constants  at  all.  Boundary  effects  between  adjacent  fiber  ribbons 
were  found  to  be  minimal,  and  determination  of  interaction  coordinates  was  robust  and 
unambiguous  because  of  the  substantial  fiber  tight  yield  (averaging  8  photoelectrons/511  keV 
interaction)  and  because  of  the  intrinsic  linearity  of  interaction  position  measurement  through  local 
tight  collection  by  fibers.  The  coordinates  obtained  were  also  accurate  to  the  very  edge  of  each 
prototype  detector,  as  opposed  to  the  Anger  coordinates  which  were  quite  distorted  near  each  edge. 

Although  our  tests  to  date  have  been  earned  out  with  R5900-L16  photomultipliers,  we  have 
also  purchased  several  R5900-M16  photomultipliers  with  16  anodes  each  4mm  x  4mm  arranged  in 
a  4  x  4  array.  These  will  be  used  to  read  out  depth-of-interaction  fibers,  and  were  intended  to  test 
readout  with  less  multiplexing  at  the  cost  of  lesser  spatial  resolution  at  the  PMT.  Very  recently, 
additional  options  have  become  available  from  Hamamatsu  Photonics  within  die  R5900  series  of 
devices;  in  particular,  they  have  produced  the  R5900-M64  photomultiplier  with  64  anodes  each 
2mm  x  2mm  arranged  in  an  8  x  8  array.  This  device  will  be  capable  of  8-fold  parallel  charge- 
division  readout  of  16  lmm-thickness  ribbons  each  with  just  2-fold  multiplexing;  another 
interesting  option  is  to  interdigitate  fibers  within  a  32mm  diameter  x  1mm  thickness  ribbon  to  make 
a  16mm  x  2mm  block,  where  the  spatial  resolution  will  still  be  dominated  by  the  upstream  optics, 
at  least  for  CsI(Na)  devices.  The  high  speed  of  LSO  permits  greater  multiplexing  without  penalty. 

7)  Data  acquisition  electronics  and  computer  interface: 

We  used  commercial  LeCroy  2249W  ADCs  to  digitize  the  outputs  of  both  4  Anger  PMTs 
per  module  and  16  multianode  PMT  charge-division  outputs  per  module.  A  trigger  was  generated 
from  the  coincidence  of  discriminated  output  of  analog-summed  Anger  PMT  signals,  and  each  of 
the  40  digitized  signals  was  delayed  by  200ns  to  land  within  a  lOOOns-width  gate.  A  second, 
lower  threshold  was  also  set  on  each  Anger  PMT  sum,  and  the  outputs  of  these  discriminators  sent 
to  a  LeCroy  2228  TDC  to  determine  the  coincidence  time  resolution  as  a  function  of  threshold. 
CAMAC  readout  of  the  42  detector  channels  was  performed  with  a  crate  controller  design  at  the 
Budker  Institute  for  Nuclear  Physics  in  Siberia,  which  is  capable  of  nearly  lOKHz  readout  of  this 
many  channels  and  is  coupled  to  a  PC.  This  rate  exceeds  that  obtainable  with  commercial  CAMAC 
interfaces  by  more  than  an  order  of  magnitude. 
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Results  and  Discussion: 


We  have  measured  the  response  of  our  prototype  mammograph  modules  and  obtained  the 
following  performance  levels  to  date: 

1 1  Spatial  resolution: 

We  have  obtained  better  than  2.5mm  FWHM  (Full  Width  at  Half  Maximum)  spatial 
resolution  in  imaging  a  point-like  source;  most  recently,  a  point  spread  function  of  less  2.0mm 
FWHM,  measured  from  a  “point  source”  which  in  fact  had  some  extension,  was  obtained 
experimentally  and  is  illustrated  in  Figure  1  on  the  following  page. 

21  Light  yield: 

We  have  obtained  an  average  of  8  photoelectrons  at  the  fiber  readout  for  1mm  fibers,  10 
photoelectrons  for  0.5mm  fibers  coupled  to  the  periphery  of  1 12mm  x  1 12mm  x  2.5mm  CsI(Na) 
crystals,  and  roughly  200  photoelectrons  at  the  Anger  PMT  readout.  The  inordinately  small  Anger 
PMT  light  yield  is  due  to  a  combination  of  several  effects:  incomplete  photocathode  coverage  by  4 
2”  round  PMTs  (we  have  purchased  and  obtained  60mm  square  PMTs  for  successor  modules), 
relatively  low  PMT  quantum  efficiency  for  the  green  wavelength-shifted  light,  and  less  than 
optimal  reflectors  about  the  module  edges  (we  plan  to  substitute  white  Tivek  reflectors  for  the 
aluminized  mylar  in  future  modules,  and  to  get  better  geometric  coverage).  We  anticipate 
increasing  the  overall  Anger  light  yield  to  roughly  400  photoelectrons/5 1 1  keV  photocapture. 

3)  Energy  resolution: 

We  have  obtained  20%  energy  resolution  FWHM,  after  correction  for  nonuniformity  in  the 
Anger  PMT  light  collection  across  the  module  (this  calibration  is  straightforward,  given  the  fiber- 
derived  event  coordinates  and  the  511keV  photopeak).  Of  this,  about  14%  FWHM  is  intrinsic 
nonuniformity,  and  about  14%  due  to  our  limited  photostatistics  (these  effects  add  in  quadrature). 
Therefore,  we  anticipate  a  future  energy  resolution  of  about  15%  FWHM  for  CsI(Na)  modules. 
This  should  be  adequate  for  clinical  scatter  rejection  requirements,  but  further  tests  are  needed. 

4)  Timing  resolution: 

We  have  obtained  50ns  FWHM  timing  resolution  for  CsI(Na)  modules  by  using  a  second 
low  threshold  on  the  Anger  sums  in  addition  to  the  threshold  used  for  energy/trigger 
discrimination.  This  is  in  agreement  with  the  decay  time  of  the  CsI(Na)  and  the  results  from 
computer  simulation  of  the  best  timing  achievable  as  a  function  of  integration  time  and  threshold 
level.  Given  the  order  of  magnitude  fast  decay  time  of  LSO,  we  would  expect  an  achieveable  LSO 
time  resolution  of  5ns  or  better.  Time  resolution  is  important  for  randoms  rejection,  particularly  at 
high  event  rates  with  efficient  detector  modules. 

5)  Rate  capability: 

Our  tests  to  date  have  been  with  inefficient  single-layer  modules,  which  have  obtained 
singles  rates  of  about  lOOKhz  and  trues  rates  of  about  lKhz  when  in  close  proximity  to  our  50 
microcurie  test  source  (10cm  separation  between  detector  modules).  The  signal-to-noise  ratio 
should  stay  about  the  same  while  our  even  rate  increases  as  we  move  on  to  higher-efficiency 
multilayer  detector  modules.  For  two  7-layer  modules  we  expect  a  nearly  50-fold  increase  in  event 
rates,  which  will  saturate  our  current  data  acquisition  (DAQ)  system  with  50  microcuries  in  the 
field  of  view.  For  this  reason,  we  are  currently  working  to  upgrade  our  DAQ  to  incoiporate 
custom  ADC  modules  with  a  PCI  interface,  designed  to  be  capable  of  up  to  100  KHZ  event  rates. 
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Figure  1  (left):  Point  Spread  Function 

Cross-section  in  x  of  the  reconstruction  of  data  from  a  disk-shaped  source  of  positrons, 
with  source  dimensions  about  6mm  diameter  by  1mm  thickness  the  thin  axis  of  the  disk  was 
oriented  along  the  x-axis.  Units  are  in  millimeters,  and  the  FWHM  (full  width  at  half  maximum)  is 
less  than  2mm.  Backgrounds  due  to  randoms  have  been  subtracted. 

Figure  2:  (right)  Tmage  of  Point  (disk)  Source 

The  source  was  oriented  as  in  Figure  1  (see  above).  The  horizontal  coordinate  is  x,  the 
vertical  coordinate  is  y,  all  units  are  in  millimeters,  and  10  contours  of  equal  intensity  are  shown. 
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Figure  3:  Image  of  Ring  Source 

In  this  image,  the  disk  source  described  in  Figure  1  and  Figure  2  was  rotated  on  a  spinning 
wheel,  with  the  disk  thin  axis  (1mm  thickness)  pointing  radially  from  the  center  of  the  circle.  All 
units  are  in  millimeters,  randoms  have  been  subtracted,  and  we  have  corrected  for  detector 
acceptance  while  accepting  all  oblique  collection  angles  which  struct  both  detectors.  The  resulting 
intensity  is  uniform  over  the  ring,  and  statistical  errors  in  background  subtraction  are  amplified  by 
the  acceptance  correction  to  large  values  only  very  near  the  edges  of  the  field  of  view. 
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Recommendations  and  Future  Plans: 


We  expect  to  build  upon  our  current  design  to  construct  prototype  detectors  with  greater 
efficiency,  higher  rate  capability,  and  improved  spatial  resolution.  Higher  efficiency  will  first  be 
obtained  by  implementing  multi-layer  CsI(Na)  modules,  starting  with  3  crystals/4  ribbons,  and 
progressing  to  7  crystals/8  ribbons.  Higher  rate  capability  will  first  be  obtained  by  substituting 
FERAs  (Fast  Encoding  and  Readout  ADCs)  for  our  current  LeCroy  2249  ADCs,  while  still 
keeping  a  CAMAC  interface  for  data  acquisition.  Later  in  the  year,  we  expect  to  fully  convert  to 
custom  ADCs  and  a  PCI  interface  which  are  currently  under  development.  Improved  spatial 
resolution  should  be  obtained  with  LSO  modules,  which  will  start  as  single-layer  modules  (the 
crystals  for  which  are  currently  being  cut  and  polished)  and  then  will  progress  to  3-layer  modules 
before  the  end  of  the  coming  year. 

Testing  of  our  modules  should  make  the  transition  from  our  Na-22  “point”  source  (which 
as  noted  earlier  is  a  6mm  diameter  disk  with  1mm  thickness)  to  a  variety  of  phantoms  which  we 
will  fill  with  F-18  is  various  liquids;  we  have  made  arrangements  to  obtain  and  transport  this 
material  from  the  cyclotron  at  the  Massachusetts  General  Hospital,  and  we  are  currently  awaiting 
radiation  safety  committee  approval  to  commence  phantom  filling  and  transportation.  The  first 
planned  “phantom”  is  hypodermic  needles  with  1mm  inner  diameter,  which  will  provide  line 
source  of  precisely  known  geometry.  The  next  phantom  will  be  a  custom-designed  object  with  a 
“DeRenzo”  pattern  (6  different  diameters  of  hot-spots  in  6  sectors  of  a  circle,  arrayed  on  the 
vertices  of  equilateral  triangles  at  separations  of  3  times  their  diameter-ranging  from  1mm 
diameters  at  3mm  separations  to  3mm  diameters  at  9mm  separations)  with  a  total  thickness  of  less 
than  5mm.  This  will  be  our  “2d”  phantom,  and  it  is  currently  in  production  at  our  machine  shop. 
Finally,  we  plan  tests  with  our  “3d”  phanto  which  is  a  commercial  “Micro  Deluxe”  phantom  we 
have  purchased  from  Data  Spectrum  corporation,  with  a  similar  pattern  but  as  cold  rods  rather  than 
hot  spots,  and  with  a  70mm  length  along  its  symmetry  axis.  Although  this  “3d”  phantom  is  in  fact 
symmetric  along  one  axis  and  thus  is  the  same  at  different  z-slices,  we  will  be  doing  3-d  imaging 
and  reconstruction  when  we  accept  events  at  oblique  incidence  and  seek  to  accurately  reconstruct 
its  entire  enclosed  volume. 

Reconstruction,  including  background  subtraction  and  perhaps  attenuation  and  scatter 
correction,  will  be  carried  out  with  a  novel  iterative  reconstruction  algorithm  which  we  have 
developed  and  which  is  described  in  Appendix  III:  “Monte-Carlo  Based  Implementation  of  the 
ML/EM  Algorithm  for  3D  PET  Reconstruction”.  Iterative  reconstruction  permits  3d  limited-angle 
tomography,  such  as  is  our  case  with  two  modules  which  do  not  entirely  surround  the  region  of 
interest  (e.g.  a  breast  or  axillary  tails).  Accurate  reconstruction  of  the  phantoms  mentioned  earlier 
are  a  necessary  prerequisite  to  later  clinical  application  of  these  or  successor  modules  (where  such 
clinical  applications  again  are  beyond  the  scope  of  this  project  as  funded  by  USAMARC). 

Finally,  we  are  working  to  cut  the  “umbilical  cord”  of  electronics  and  mechanical  mounts 
which  currently  restrict  our  test  modules  to  our  laboratory,  and  to  construct  new  modules 
connected  to  portable  electronics  and  a  portable  data  acquisition  system,  mounting  on  our  donated 
mammography  stand.  Our  intention  is  to  transport  our  prototype  modules  to  a  clinical  setting  at  the 
earliest  possible  date,  and  if  possible  we  will  perform  some  simple  measurements  in  a  clinical 
setting  to  ascertain  the  radiation  environment  which  these  modules  will  encounter  in  clinical 
practice.  We  anticipate  simulating  such  environment  by  positioning  distributed  activity  throughout 
the  regions  outside  the  field  of  view  during  phantoms  studies  toward  the  end  of  next  year,  as  we 
approach  the  end  of  this  program. 
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Conclusions 


Because  of  the  limitations  of  standard  X-ray  mammography  for  the  early  detection 
of  breast  cancer,  in  particular  for  women  with  radiographically  dense  breasts,  nuclear  medicine 
imaging  techniques  are  of  interest  for  this  application.  Of  these  techniques,  Positron-Emission 
Tomography  has  the  greatest  potential  for  high  sensitivity  to  small  tumors,  and  good  uptake  ratios 
have  been  demonstrated  for  breast  tumors  with  the  standard  PET  radiotracer  FDG.  The  key  PET 
detector  attributes  for  this  application  are  high  acceptance  and  good  efficiency  at  high  rates,  fine 
spatial  resolution,  and  good  scatter  rejection.  We  are  constructing  and  testing  novel  prototype 
PET  detector  module  which  have  demonstrated  high  spatial  resolution  (2.5mm  FWHM  in  all  three 
dimensions  when  reconstructing  a  point  source),  with  large  acceptance  and  a  unique  active-edge 
and  active-comer  geometry  which  will  be  compatible  with  standard  breast  imaging  procedures 
including  the  imaging  of  axillary  nodes.  We  plan  to  continue  to  optimize  and  test  larger-scale, 
higher-efficiency,  higher-speed  and  higher-resolution  modules  within  the  coming  year. 
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Abstract 

By  using  wavelength-shifting  libers  coupled  to  thin  plates 
of  LSO  and  to  photomultipliers,  we  have  demonstrated  spatial 
resolution  of  2mm  FWHM  tor  photocapture  events. 
Optimized  systems  are  believed  to  be  capable  ol  better  than 
1mm  FWHM  resolution,  without  requiring  either  very  small 
crystals  or  large  numbers  of  photosensors.  With  this  sensor 
resolution,  annihilation  acoilinearity  dominates  system 
resolution  for  large  PET  rings.  By  radially  stacking 
alternating  layers  of  thin  crystals  and  wavelength-shifting 
fibers,  a  direct  depth-of-interaction  measurement  is  possible. 
Energy  measurements  and  coincidence  timing  may  be 
provided  by  trigger  photomultipliers  in  a  dual-photodetector 
geometrv.  Results  from  computer  simulations  and  component 
response  measurements  on  proot-ol-conccpt  prototypes  are 
presented. 

I.  INTRODUCTION 

We  are  developing  very  high-precision  PET  detector 
modules  with  depth-of-interaction  sensitivity,  using  a  new 
design  incorporating  wavelength-shifting  liber  readout  ol  thin 
LSO  crystals.  By  building  modules  from  thin  crystal  layers 
stacked  radially,  one  can  obtain  both  depth-ot-interaction 
information  and  fine  spatial  resolution  in  the  tangential  and 
axial  directions.  Component-level  testing  of  this  technique 
has  demonstrated  <2  mm  FWHM  spatial  resolution  lor  2mm 
thick  LSO  crystals.  Optimization  and  extrapolation  to  full 
PET  systems  may  yield  devices  with  resolution  as  small  as 
1mm  FWHM  using  l-2mm  thick  crystal  layers:  the  resolution 
is  expected  to  vary  linearly  with  the  crystal  layer  thickness. 
Depth-of-interaction  sensitivity  is  essential  for  uniformly 
high-precision  PET  throughout  a  large  tield-ot-view.  it  one  is 
to  preserve  the  detector  thickness  needed  for  high  efficiency 
and  high  event  statistics.  Depth-of-interaction  sensitivity  also 
minimizes  the  necessary  detector  diameter,  thereby  decreasing 
both  acoilinearity  errors  and  system  costs,  while  increasing 
solid  angle  coverage.  Our  modular  device  design  permits 
testing  of  a  variety  of  detection  geometries,  including 
geometries  appropriate  tor  laboratory  animal,  human  brain, 
and  dedicated  PET  mammography  applications. 

II.  MATERIALS  AND  METHODS 

We  have  used  ribbons  of  wavelength-shifting  fibers  to 
determine  the  positions  of  gamma-ray  interactions  within  thin 
LSO  scintillator  crystals.  Within  a  wvaelength-shifting  fiber, 
fluorescent  dopant  molecules  absorb  short-wavelength 
incident  photons  and  then  isotropically  emit  secondary  longer- 
wavelength  photons.  One  can  use  the  fluorescent  re-emission 


process  effectively  to  inject  light  into  a  light  pipe  through  the 
sides  of  the  pipe,  since  a  fraction  of  the  re-emitted  light  is 
trapped  and  piped  along  the  fiber.  In  this  way,  the  fiber 
becomes  a  light  sensor,  with  short-wavelength  light  as  input 
and  long-wavelength  light  as  output.  An  array  of  parallel 
fibers  (i.e.  a  fiber  ribbon)  may  then  be  used  to  measure 
incident  light  profiles  in  one  dimension;  two  perpendicular 
ribbons  (on  opposite  faces  of  the  thinscintillator  crystal)  can 
provide  a  two-dimensional  incident  light  profile  measurement. 
Gamma-ray  interactions  in  dense  scintillators  result  in  a  very 
nearly  point-like  light  source,  because  of  the  short  range  of 
secondary  photocapture  electrons.  For  polished  high-index 
scintillator  crystals,  only  light  emerging  nearly  normal  to  the 
crystal  surface  escapes  total  internal  reflection,  resulting  in  an 
emission  spot  with  diameter  comparable  to  the  crystal 
thickness.  We  earlier  described  this  technique  as  applied  to 
BGO  crystals  read  out  with  red  wavelength-shifting  fibers  [1]. 

In  the  present  work  we  have  focussed  on  the  fast,  high-Z, 
high-brightness  scintillator  Lutetium  Oxyorthosilicate  or  LSO 
[2].  We  have  performed  most  of  our  initial  measurements  on 
two  samples  of  LSO,  each  2mm  x  10mm  x  1 0mm,  which  we 
obtained  from  Schlumberger-Doll  Research  [3].  These 
crystals  were  hand-polished  on  their  large  surfaces,  and  all 
surfaces  not  coupled  to  fibers  were  painted  black.  Fibers  were 
coupled  to  crystals  with  a  high-index  epoxy  (n=l  .65),  EpoTek 
OG127  [4],  For  readout  fibers  we  used  1-mm  diameter, 
round,  doubly-clad  wavelength-shifting  fibers  from  Kuraray 
Corp  [5];  we  obtained  our  best  spectral  matching  with  their 
spectral  type  Y- 1 1 .  Doubly-clad  Fibers  have  a  core  refractive 
index  of  1 .60,  an  inner  cladding  index  of  1 .49,  and  an  outer 
cladding  index  of  1.42;  this  yields  a  numerical  aperture  of  0.74 
and  a  trapping  efficiency  of  about  7 %  per  fiber  end  for  the  re¬ 
emitted  light  (the  remaining  86%  of  the  re-emitted  light  is  not 
trapped  and  emerges  through  the  sides  of  the  fiber  ribbon). 

Our  initial  tests  have  been  carried  out  with  the  apparatus 
illustrated  schematically  in  Figure  1.  We  form  a  pencil  beam 
of  51 1-keV  gamma  rays  by  requiring  the  coincident  detection 
of  two  back-to-back  annihilation  gammas  produced  by  a 
nearly  point-like  Na-22  positron  source.  The  first  of  the  two 
gammas  strikes  a  small  1mm  x  3mm  x  15mm  LSO 
“triggering”  crystal  through  its  1mm  x  3mm  entry  face,  while 
the  partner  ray  goes  through  a  1mm  lead  collimator  to  the 
“target”  2mm  x  10mm  x  10mm  crystal.  The  “target”  crystal  is 
coupled  to  a  ribbon  of  10  parallel  1mm  diameter  round 
wavelength-shifting  fibers  and  to  a  trigger  PMT.  Since  the 
light  collected  by  the  trigger  PMT  is  proportional  to  the 
energy  deposited  in  the  target  crystal,  we  are  able  to  select 
events  with  photocapture  of  a  pencil-beam  gamma  by  the 
target  crystal,  by  requiring  a  large  trigger  signal  in  coincidence 
with  a  signal  from  the  “triggering”  crystal.  The  "triggering” 
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LSO  crystal  and  its  PMT.  the  Na22  source  and  collimator  are 
all  mounted  on  a  calibrated  X-Y  adjustable  table,  while  the 
“target”  crystal  and  PMT  with  fiber  assembly  are  fixed:  in  this 
way.  a  pencil  beam  can  be  precisely  directed  into  the  target 
setup,  and  may  be  swept  along  the  target.  All  measurements 
are  made  with  unshaped  signals  directly  coupled  to  charge- 
sensitive  ADCs  (LeCroy  2249W);  the  timing  gate  width  was 
150nsec  and  only  photopeak  events  were  used  for  analysis. 


Figure  I :  Schematic  illustration  of  test  apparatus. 

Each  of  10  fibers  coupled  to  the  target  crystal  is  read  out 
by  an  individual  Hamamatsu  R1635  PMT.  allowing  us  to 
record  the  light  profile  and  centroid  position  within  the  fiber 
ribbon  on  an  event-bv-event  basis:  all  PMTs  were  calibrated^ 
and  outputs  converted  to  photoelectrons.  The  index  of 
refraction  of  LSO  is  about  1.8.  so  we  used  the  high-index 
optical  epoxy  for  coupling  crystals  to  libers.  The  thin  sides 
were  blackened  with  an  opaque  mixture  of  high-index  epoxy 
and  xerographic  toner  (to  simulate  a  larger  crystal).  Fibers 
were  coupled  to  both  of  the  large  LSO  surfaces,  with  30cm 
long  fibers  having  one  polished  end  at  the  readout  PMTs  and 
their  other  ends  1  cm  from  the  crystal.  Light  trapped  at  the 
optical  interlace  between. cladding  and  air  (  cladding  light  ) 
was  removed  by  an  opaque  fiber  holder  gripping  each  fiber 
before  its  coupling  to  its  readout  PMT.  The  tiber/PMT  optical 
interface  was  made  with  optical  grease. 

We  have  developed  a  complete  Monte  Carlo  simulation 
of  our  detector  s  optics,  including  all  optical  interfaces  and  the 
correct  geometry  for  round  multiclad  libers.  We  have  also 
developed  a  detailed  radiation  transport  Monte  Carlo 
simulation  to  study  the  effects  of  gamma-rav  scattering  and 
attenuation  both  in  objects  being  imaged  and  within  the 
detector  itself.  The  latter  simulation  is  based  on  the  GEANT 
detector  simulation  code  commonly  used  in  high-energy 


physics,  which  allows  for  simple  implementation  and 
modifications  of  the  detector  geometry  through  a  flexible  and 
powerful  user  interface.  We  generate  gamma-rays  from  our 
source  phantom  according  to  a  positron  range  distribution 
appropriate  to  18F,  and  with  acollinearity  between  the 
annihilation  gammas  sampled  from  a  Gaussian  distribution 
with  width  0.3  degrees.  The  detector  resolution  was  usually 
taken  to  match  the  (not-yet-optimized)  response  functions 
measured  with  our  test  apparatus,  and  matches  the  measured 
dependence  on  light  yield  (gamma  ray  energy). 

We  have  also  begun  work  on  making  images  with 
simulated  data  using  a  simple  2-D  filtered  back-projection 
reconstruction  algorithm,  for  a  number  of  test  phantoms.  In 
particular,  we  have  looked  at  the  benefits  of  depth-of- 
interaction  resolution  for  small-diameter  systems  such  as  for 
use  with  laboratory  animals.  For  some  of  these  studies  we 
have  simulated  a  1 /2-scale  version  of  a  standard  resolution 
phantom,  with  resulting  10cm  diameter  axially  symmetric 
phantom  viewed  by  a  simulated  20cm  diameter  detector  ring. 
To  contrast  with  our  detector,  we  have  modelled  a  device 
without  depth-of-interaction  sensitivity,  consisting  of  8mm 
thick  LSO  on  a  20cm  diameter  ring.  For  the  radially- 
segmented  fiber  readout  20cm  diameter  ring  we  have  used  a 
total  crystal  thickness  of  24mm,  which  would  result  in 
improve  coincident  efficiency  by  a  factor  of  3.2  relative  to  a 
8mm  thickness  ring. 

A  wide  range  of  system  design  options  are  open  with  a 
radially-segmented  fiber  readout  geometry,  and  we  have 
begun  our  exploration  of  these  options  with  a  baseline  module 
having  thickness  24mm  and  area  128mm  x  128mm.  Our  plan 
is  to  construct  and  test  at  least  one  such  module  within  the 
coming  year.  Assuming  2mm  thick  LSO  layers  and  1mm 
diameter  fibers,  this  requires  the  readout  of  13  x  128  =  1664 
fibers  per  module.  Multiplexing  is  clearly  necessary  in  order 
minimize  the  cost  for  readout  photomultipliers  and  digitizing 
electronics.  One  method  for  such  multiplexing,  employing 
separate  readout  of  two  fiber  ends,  is  illustrated  schematically 
in  Figure  2  below.  Other  schemes,  including  mirrored  fibers 
and  two  different  colors  of  waveshifting  fiber,  are  possible. 


One  Plane  of  8  Parallel  Fiber  Ribbons 


Profile  across  each  of  Index  for  each  of 

8  fiber  ribbons  8  fiber  nbbons 

multiplexed  onto  16  anodes  from  8  anodes 


Figure  2:  Multiplexed  readout  of  fiber  Ribbons  using  16- 
channel  multianode  photomultipliers. 
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We  have  been  testing  a  new  R5800-L16  multianode 
photomultiplier  from  Hamamatsu  Corporation,  which 
combines  16  input  pixels  each  1mm  x  16mm  in  area  with  a 
very  narrow  pulses  (2ns  FWHM)  and  good  single¬ 
photoelectron  resolution  [6].  By  reading  out  the  two  ends  ot 
multiclad  fiber  ribbons  to  two  different  multianode  PMTs. 
one  can  read  out  as  many  as  256  liber  ends  per  PMT, 
Multiplexed  readout  of  fiber  ribbons  corresponding  to 
different  crystal  layers  is  also  possible  and  probably  desirable. 
If  0.5mm-diameter  fibers  are  used,  the  physical  limit  is  1000 
fiber  ends  per  PMT.  Although  not  indicated  in  the  figure,  it 
may  be  desirable  to  break  up  each  128mm  x  128mm  crystal 
plane  into  64  16mm  x  16mm  optically  isolated  >.ells  to 
reduce  multiplexing  cross-talk. 

Electronics  channel  count  can  be  minimized  by  charge- 
division  readout  (two  signals  rather  than  16  trom  each 
photomultiplier)  as  is  typically  done  with  crossed-wire 
position  sensitive  PMTs.  We  have  measured  20ns  pulse 
widths  (into  a  50  ohm  load)  and  «lmm  resolution  at  20 
photoelectrons  for  charge  division  readout  ot  our  R5800-L16 
PMT.  using  100  ohm  resistors  to  connect  the  anodes.  In 
principle,  the  readout  ol  a  16-laver  module  with  area  1-8  mm 
x  128mm  using  0.5mm  libers  could  require  tust  4  multianode 
photomultipliers  (at  a  current  cost  of  SI 500  each),  with  8 
associated  fast  ADC  channels. 


III.  PRELIMINARY  RESULTS 

Our  measured  light  yield  for  fiber  readout  of  LSO  crystals 
is  in  general  agreement  with  our  expectations.  The  light 
output  for  LSO  is  quoted  at  venous  levels  in  (he  literature,  and 
is  a  strong  function  of  the  concentration  of  Cerium  dopant:  the 
most  recent  work  cites  a  light  yield  ot  20.000  photons/MeV 
deposited  or  about  10.000  photons  per  annihilation  gamma 
photocapture  event  tor  LSO  with  a  Cerium  concentration  ot 
0.227r  (7).  Monte  Carlo  simulations  predict  that  more  than 
259!-  of  emitted  photons  should  emerge  from  a  given  crystal 
readout  face  into  the  corresponding  liber  ribbon,  with  photons 
incident  at  large  angles  (relative  to  the  surface  normal)  totally 
internally  reflected.  Direct  coupling  of  a  photomultiplier  with 
20 7r  quantum  efficiency  (when  averaged  over  the  LSO 
spectrum)  to  LSO  with  this  brightness  should  thus  give  about 
500  photoelectrons:  this  is  consistent  with  our  measured  light 
vield  in  this  "direct  coupled”  geometry.  The  next  link  in  the 
optical  chain  is  the  match  of  the  emission  spectrum  ot  the 
LSO  crystal  to  the  absorption  spectrum  of  the  fiber.  This 
factor  is  about  609!-.  based  on  convolution  of  published  LSO 
emission  spectra  and  Kurarav's  absorption  spectra  for  Y-l  1. 
We  thus  expect  that  roughly  1500  scintillation  photons  are 
absorbed  and  re-emitted  within  the  wavelength  shifting  fibers 
in  each  photocapture  event. 

Multiclad  fibers  trap  149!-  of  their  re-etnitlcd  light,  or  77c 
per  fiber  end.  Our  measurements  comparing  light  yield  with 
multiclad  vs.  singly  clad  fibers,  where  the  latter  have  47c 
trapping  per  end.  confirms  this  trapping  fraction  tor  the 
multiclad  fibers  we  used.  This  implies  that  with  multiclad 
fibers  about  100  photons  are  piped  toward  each  fiber  end  each 
photocapture  event.  Finally,  some  photons  were  attenuated 
within  the  readout  fibers  (such  self-absorption  filters  out  the 


shortest- wave  length  photons,  and  typically  results  in  re¬ 
emission  of  an  unpiped  green  photon),  and  a  few  photons 
arriving  at  the  photomultiplier  were  reflected  at  the  optical 
interface  between  fiber  and  photomultiplier.  Our  estimates  of 
these  effects  were  that  we  had  about  80  photons  entering  our 
fiber  readout  photomultipliers  each  photocapture  event,  where 
they  had  an  average  quantum  efficiency  of  about  107c  for 
being  converted  to  photoelectrons  (Y-l  1  has  its  emission  peak 
at  495  nm).  This  accounts  for  the  roughly  8  photoelectrons 
which  we  observed  per  photocapture  event  with  our  test 
apparatus.  This  light  yield,  although  small,  is  sufficient  for 
multiplexed  readout  with  good  efficiency.  Additional  light 
yield  may  be  achievable  by  using  fibers  with  a  better  spectral 
match  to  LSO  emission,  by  using  PMTs  with  better  quantum 
efficiency  for  green  light  (green-extended  R5800-L16  PMTs 
are  currently  under  development  at  Hamamatsu)  or  by 
obtaining  brighter  LSO. 

The  emission  light  profile  from  LSO  crystals  and  the 
measured  spatial  resolution  achieved  at  our  current  light  level 
is  shown  in  Figure  3  below.  The  resolution  horizontal  scale 
has  been  calibrated  by  moving  the  pencil  beam  in  our  test 
fixture  across  the  crystal  by  known  increments  and  measuring 
the  resulting  changes  in  the  measured  centroid  position. 


Figure  3:  Measured  response  of  LSO/Fiber  test  apparatus. 

Left  =  Predicted  and  measured  profile  of  light  collection. 

Right=  Measured  spatial  resolution  using  centroids. 

By  measuring  the  distribution  of  photoelectrons  across 
the  fiber  ribbon  for  photocapture  events,  we  were  able  to 
reconstruct  the  event  position  with  an  accuracy  of  about 
0.7mm  RMS.  There  is  still  a  contribution  from  the  beam 
width  to  this  resolution,  for  which  we  have  not  corrected.  The 
accuracy  of  our  coordinate  measurement  varies  directly  with 
the  width  of  the  light  emission  profile,  so  that  the  quality  ot 
crystal  surface  finish  is  important.  Our  Monte  Carlo  optical 
simulation  predicted  a  profile  of  about  1.0mm  RMS  for  a 
2mm  thick  ideally  polished  LSO  crystal,  while  the  measured 
profile  had  a  width  of  1.6mm  RMS.  Evidently,  our  crystal 
surface  preparation  and  fiber  coupling  could  be  improved. 
Spatial  resolution  improves  with  light  yield  like  the  square 
root  of  the  number  of  photoelectrons  N,  and  ideally  the 
resolution  should  approach  the  profile  width  /  'i  N. 
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We  have  carried  out  a  simulation  tor  a  geometry 
appropriate  to  Positron-Emission  Mammography,  which  we 
have  simplified  as  two  parallel  planes  separated  by  I Ocm  of 
water.  Figure  4  illustrates  the  line-spread  function  (LSF) 
predicted  for  a  perfect  detector  (due  to  just  positron  range  and 
acollinearity),  for  fiber  readout  from  1mm  thick  LSO.  and  for 
fiber  readout  from  2mm  thick  LSO:  the  latter  two  are 
extrapolated  from  our  non-optimized  test  results.  It  is 
apparent  that  range  and  acollinearity  effects  are  very  small . 


As  mentioned  earlier,  we  have  contrasted  the  simulated 
response  of  our  proposed  detector  with  that  ot  a  device  with 
comparable  spatial  resolution  but  without  depth-ot-interactiori 
sensitivity;  both  were  formed  into  20cm  diameter  rings. 


Figure  5  shows  the  reconstructed  resolution  as  a  function 
of  distance  from  the  center  of  the  field-of-view  for  each 
device).  Again,  we  have  used  the  measured  but  not-yet- 
optimized  resolution  from  our  test  apparatus  as  input  to  this 
simulation;  in  addition,  we  have  indicated  the  resolution 
which  should  be  achievable  with  1mm  thick  LSO  layers  with 


our  fiber  readout.  While  our  device  has  very  nearly  uniform 
resolution  over  the  entire  field  of  view,  the  device  without 
depth-of-interaction  sensitivity  has  clear  resolution 
degradation  away  from  the  center  of  the  field  of  view. 

IV.  DISCUSSION 

Radially  segmented  readout  of  thin  LSO  crystals  with 
wavelength-shifting  fibers  has  the  following  advantages: 

•  Very  High  Spatial  Resolution 

•  Depth-of-Interaction  Measurement 

•  High  Rate  Capability 

•  Potential  for  Low  Cost 

•  Flexible  Readout  Options 

Unlike  most  designs  for  very  high-resolution  PET  detectors, 
this  design  scales  readily  to  large  systems  while  preserving 
high  efficiency  and  fine  resolution;  it  also  only  makes  use  of 
existing  photosensor  and  electronics  technology. 

Since  the  spatial  resolution  scales  linearly  with  the  crystal 
thickness,  extremely  fine  spatial  resolution  should  be 
achievable  with  small-diameter  rings.  Resolution  blurring  due 
to  acollinearity  (about  2.2mm  FWHM  per  meter  of  ring 
diameter)  can  be  made  to  dominate  system  resolution  for  large 
rings;  smaller  rings  can  achieve  submillimeter  resolution  with 
high  efficiency  across  a  large  field  of  view  (if  short-range 
positrons  like  18F  are  used)  by  using  very  thin  crystal  layers. 
Axial  resolution  is  naturally  comparable  to  tangential 
resolution,  and  to  depth-of-interaction  resolution.  Double-hits 
from  Compton-scatters  within  the  detector  (which  can  degrade 
resolution  of  block  detectors)  are  readily  identifiable  with 
suitable  electronics,  since  they  will  almost  always  deposit 
energy  in  more  than  one  crystal  layer.  Use  may  be  made  of 
Compton  scattered  (rather  than  photocaptured)  events  in  small 
rings  imaging  objects  with  little  scatter,  with  only  a  small 
degradation  in  spatial  resolution.  Finally,  detector  modules 
may  be  combined  in  a  “gapless"  geometry  by  radially 
overlapping  in  a  pinwheel  arrangement  as  shown  in  Figure  5, 
any  oversampling  which  which  results  in  the  overlap  regions 
may  be  corrected  using  depth-of-interaction  information. 


200  mm 


Figure  5:  Conceptual  System  Geometries  with  Fiber  Readout 
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Depth-of-interaction  measurement  is  critical  for 
preserving  very  fine  spatial  resolution  throughout  a  large  field 
of  view,  without  compromising  on  detector  coincidence 
efficiency  (which  is  critical  for  the  statistics  needed  to  made  a 
high-resolution  image).  Depth-of-interaction  sensitivity 
permits  a  reduction  in  ring  diameter,  which  pays  several 
dividends:  lower  system  costs  for  less  scintillator  and  fewer 
readout  channels,  smaller  effects  from  acollinearity,  and 
higher  efficiency  through  greater  solid  angle  coverage.  The 
last  comes  at  the  price  of  increased  scatter  acceptance,  but  this 
trade-off  is  already  being  chosen  in  the  selection  of  3-D  over 
2-D  PET.  Increasing  the  module  size  (both  in  axial  and 
tangential  dimension)  yields  a  more  cost-effective  readout  by 
using  longer  readout  fibers,  which  increases  the  detector 
volume/PMT  area  ratio;  the  price  paid  is  lower  rate  capability. 

With  fully  multiplexed  readout  (including  in  depth)  of 
128mm  x  128mm  area  modules,  a  rate  capability  of 
>lMhz/moduie  should  be  achievable.  The  45ns  decay  time  of 
LSO  permits  signal  integration  over  150ns,  with  digization  of 
charge-division  signals  to  the  required  7-bit  resolution 
possible  in  less  than  1  microsecond.  Time  correspondence  of 
8  PE  signals  from  45ns  scintillator  has  high  efficiency  within 
a  10-20ns  window.  The  coincidence  timing  (for  singles 
rejection)  and  energy  resolution  achievable  with  an  all-fiber 
readout  are  more  problematic.  If  necessary,  fiber  readout  can 
be  supplemented  with  large  trigger  PMTs  which  would  collect 
the  re-emitted  but  non-trapped  light  which  emerges  through 
the  sides  of  the  wavelength-shifting  fibers  and  to  which  both 
fibers  and  crystals  are  transparent.  Our  measurements  have 
shown  that  although  less  than  \0%  of  LSO  light  passing 
through  one  fiber  layer  is  absorbed  in  the  next,  the  layers 
which  are  opaque  to  LSO  light  are  very  transparent  to  fiber 
light.  Supplemental  readout  with  large  trigger  PMTs  would 
yield  hundreds  of  photoelectrons,  sufficient  for  Ins 
coincidences  and  energy  resolution  limited  by  etfects  other 
than  photostatistics.  Large-area  photodiodes  might  be 
substituted  for  photomultipliers  to  preserve  compatabiiity  with 
operation  in  a  magnetic  field,  if  necessary. 

Fiber  readout  of  LSO  is  potentially  quite  cost-effective, 
presuming  that  the  bulk  price  of  LSO  continues  to  fall.  It  is 
currently  at  $50/cc  in  quantities  of  100’s  of  cc's,  and  is 
projected  by  the  crystal  grower  to  fall  to  S15/cc  in  about  1 
year.  At  this  price  it  becomes  competitive  with  BGO,  and 
fiber  cutting  and  polishing  costs  become  signficant  (especially 
for  very  high-resolution  devices  with  small  crystals).  The 
crystal  cutting  and  polishing  required  for  the  thin  Hat  crystals 
used  with  fiber  readout  is  significantly  less  than  for  the  long 
thin  crystals  in  other  detectors,  and  we  will  carefully  check 
how  our  performance  varies  with  crystal  surface  treatment. 
The  cost  for  4  multianode  PMTs  to  read  out  fibers  from  a 
128mm  x  128mm  module  is  comparable  to  the  25  1"  PMTs 
which  would  be  required  to  read  out  the  same  area  in  a  block 
detector;  if  larger  modules  have  sufficient  rate  capability, 
fiber  readout  costs  fall  still  further.  The  cost  of  fibers  in 
quantity,  about  $  1/in  for  0.5mm  diameter,  makes  the  fiber  cost 
small  in  comparison  to  crystals  and  PMTs.  Finally,  the 
electronics  channel  count  reduction  yields  savings,  reducing 
the  25  ADCs  for  the  above  block  detector  to  8  in  our  design. 


Although  our  initial  results  have  been  encouraging,  we 
have  been  somewhat  handicapped  by  the  availability  of  only 
small  samples  of  LSO.  This  situation  has  now  changed,  and 
we  anticipate  the  availability  of  hundreds  of  cc's  of  LSO  by 
early  1996.  Our  plan  is  to  further  optimize  and  to  construct 
large  modules  with  this  new  material  as  soon  as  possible. 

V.  CONCLUSIONS 

We  have  demonstrated  that  wavelength-shifting  fibers 
can  be  used  to  read  out  gamma-ray  interaction  positions 
within  thin  sheets  of  LSO  crystals,  with  high  efficiency.  The 
spatial  resolution  achievable  (FWHM)  is  less  than  the 
thickness  of  the  crystal  layer  being  read  out.  Consequently, 
by  constructing  a  detector  from  a  thick  stack  of  many  thin 
layers,  one  can  combine  very  fine  spatial  resolution  with 
depth  of  interaction  sensitivity  and  high  efficiency.  With 
optical  multiplexing  at  the  crystal  layers  and  at  the  readout 
photomultipliers,  a  cost-effective  design  with  extremely  high 
spatial  resolution  is  achievable;  this  design  is  readily  scaled  to 
large-diameter  systems.  During  the  coming  year,  we  plan  to 
constract  one  or  more  128mm  x  128mm  x  20mm  thick 
modules  to  test  and  optimize  this  design. 
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Abstract 

We  have  characterized  the  new  Hamamatsu  R5900-L16 
multianode  photomultiplier  with  regard  to  its  suitability  for 
position-sensitive  readout  with  just  two  analog  channels  for 
each  16-anode  PMT.  Because  of  the  low  capacitance  of  the 
anode  structure  of  the  PMT,  we  obtain  an  output  pulse  width 
of  less  than  20  ns  for  narrow  input  pulses.  W’e  were  able  to 
obtain  position-sensitive  readout  of  single  photoelectron 
pulses  with  a  spatial  width  of  less  than  l  mm  FWHM  for  a 
collimated  light  source  and  about  2  mm  FWHM  for  single 
fiber.  Combined  with  the  good  single  phoioclcctron  peak  for 
R5900-L16,  this  provides  excellent  performance  at  very  low 
light  levels.  We  discuss  application  of  this  device  to  the 
precise  determination  of  gamma-ray  interaction  position 
within  thin  crystal  scintillators. 

I.  Introduction 

Position  sensitive  photomultipliers  have  been  used  by  a 
number  of  groups  for  the  readout  of  scintillating  or 
wavelength-shifting  fibers.  Crossed-wirc  and  multianode 
PMTs  have  also  been  used  for  the  position-sensitive  readout  of 
extended  light  sources  [2-5].  At  low  light  levels,  the  cross¬ 
talks  within  the  PMT  can  result  in  significant  degradation  in 
spatial  resolution.  There  are  several  sources  ot  such  cross¬ 
talks:  light  spreading  through  the  entry  window  glass  on  the 
way  to  the  photocathode,  imperfect  electron  optics  between 
photocathode  and  dynode  structure,  secondary  electron  shower 
spreading  within  the  multiplier  assembly,  or  just  electrical 
cross-talks  between  anodes  (induced  signals).  A  significant 
difference  between  position-sensitive  cross-wired  PMTs  like 
the  Hamamatsu  R2486,  versus  multianode  PMTs  like  the 
Hamamatsu  R4760,  is  the  considerably  higher  cross-talks  and 
therefore  decreased  spatial  resolution  at  low  light  level  of  the 
former.  The  Hamamatsu  R2486  for  instance  requires  10000 
photons  at  490  nm  (about  1000  photoelectrons)  to  obtain  a 
spatial  resolution  of  0.2  mm  FWHM  [2],  the  resolution  is 
proportional  to  the  square  root  of  the  number  of  photoelectron 
per  event.  Another  difficulty  with  cross-wired  PMTs  for  some 
applications  is  its  limited  rate  capability  due  to  significant 
pulse  width  with  typical  resistor  chain  readout. 
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Nonetheless,  the  reduced  electronics  channel  count  and 
lower  device  cost  (particularly  per  unit  of  useful  photocathode 
area)  associated  with  crossed-wire  PMTs  have  made  them 
attractive  for  several  applications,  including  for  the  readout  of 
sets  of  scintillating  fibers. 

There  is  some  success  in  attempts  to  use  charge  division 
method  for  two-dimensional  readout  of  an  array  of  optical 
fibers  by  multianode  PMT  having  just  four  analog  channels 
[5].  But  again,  one  needs  hundreds  of  photons  to  obtain 
submillimeter  resolution. 

We  investigated  the  properties  of  the  newly  developed 
PMT  -  the  Hamamatsu  R5900-L16  -  because  of  its  very 
promising  specifications  caused  mainly  by  recently  developed 
microchannel  dynode  structure. 

II.  Experimental  Setup 

A.  R5900-L16  -  multianode  PMT 

The  photomultiplier  described  here  has  external 
dimensions  of  merely  a  cube  25*25*25  mm.  It  has  flat 
bialkali  photocathode,  10-stage  dynode  structure  and  16  anode 
stripes  0.8  *  16  mm  regularly  spaced  at  1  mm  pitch.  An 
important  feature  of  the  PMT  distinguishing  it  from  previous 
ones  is  that  the  microchannel  dynodes  work  effectively  like 
separate  PMTs  in  terms  of  electron  optics.  This  minimizes 
electron  shower  smearing  and  allows  one  to  obtain  very  fine 
spatial  resolution  as  will  be  shown  below.  Another  useful 
feature  of  the  microchannel  structure  is  a  relatively  high  gain 
at  low  operating  voltage  -  about  3..4*106  at  800V  .  The  first 
prototype  of  R5900-L16  sent  to  us  by  Hamamatsu  had  a  size 
of  photocathode  of  15.5  *  15.5  mm  which  is  a  0.5mm 
narrower  than  dynode  structure.  Most  of  the  measurements 
were  done  on  this  particular  PMT.  The  photocathode  size  of 
the  last  version  of  R5900-L16  is  significantly  larger  - 
approximately  21  *  22  mm.  We  have  repeated  some  of  the 
measurements  to  check  the  impact  of  this  design  change  on  the 
performance. 

B.  Charge  Division  Readout 

We  have  connected  16  anodes  of  the  PMT  by  a  linear  array 
of  15  100-ohm  resistors  to  couple  16  individual  outputs  into  2 
charge  division  outputs.  Each  of  the  two  signals  was  directly 
coupled  through  50-ohm  cables  to  a  LeCroy  2249W  operating 
with  a  150  ns  gate  width.  To  minimize  artifacts  caused  by 
discrete  nature  of  ADC  data  we  have  inserted  10-fold 
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gnplificadon  of  the  charge  division  signals  giving  us  a  single 
piotoelectron  peak  at  about  35  channels  of  the  ADC  over  the 
fcdestal. 


I 


fj.  Data  Acquisition  System 

1  For  measurements  described  here  we  used  rather  simple 
data  acquisition  system  consisting  of  a  CAMAC  ADC 
mentioned  above  and  micro-VAX  computer  under  Open  VMS 
ponnected  to  the  CAMAC  crate  and  controller  via  SCSI-bus. 
Data  collected  on  a  disk  of  the  computer  then  were  analyzed  in 
it  PAW  (CERN  written  data  analysis  software)  session. 


i  III.  Tests  and  Measurements 

£ 

%  Timing  and  Single  Photoelectron 
?  Characteristics 

The  PMT  itself  (when  its  anodes  are  read  out  separately) 
has  fast  rise  time  and  narrow  pulse  width  -  0.6  ns  nse  time.  1 .7 
ns  fall  time  when  loaded  to  50  Ohm.  After  connecting  all 
anodes  by  a  linear  array  of  100-Ohm  resistors,  we  obtain  a 
pulse  width  of  less  than  20  ns.  To  check  single  photoelectron 
amplitude  peak  we  connected  together  two  ends  of  the  charge 
division  chain.  The  PMT  was  illuminated  by  a  green  LED 
giving  0.1  photoelectron  in  average.  The  distribution  is  shown 
on  Fig.  1.  The  peak-to- valley  ratio  of  almost  3:1  is  a  very 
desirable  characteristic  for  low  light  level  applications. 


Amplitude,  cnc°°e<  ADC 


Figure  1:  Single  photoelectron  peak  of  the  R5900-L16.  All  anodes 
are  read  out  together.  HV=800V. 

B.  Charge  Division  Calibration 

Ideally,  if  we  know  parameters  of  the  chain  -  we  can 
calculate  the  calibration  coefficient  between  charge  division 
ratio  R=(A2-A1)/(A1+A1)  (A1  and  A2  -  amplitudes  measured 
from  two  ends  of  the  chain)  and  position  in  the  chain, 
considering  each  anode  to  be  a  current  source.  The 
correspondence  between  charge  division  node  number  and 
measured  ratio  R  should  be  linear.  Therefore  even  if  several 
photoelectrons  are  distributed  over  several  PMT  anodes  the 
measured  ratio  R  would  correspond  to  the  center-of-gravity  of 


received  photoelectrons.  In  our  particular  case  we  have 
Rin=50  Ohm  and  15  resistor  of  the  chain  100  Ohm  each.  So 
ideally  for  point  source  from  the  first  or  the  last  anode  we 
should  measure  the  ratio 

R=+-(50+ 15*1 00-50)/(50+ 15*1 00+50)=+-0.9375 . 

In  real  situation  one  has  to  correct  this  ratio  for  different 
gains  and  input  impedances  of  the  two  amplifiers.  Practically 
we  found  that  the  ratio  R  ranges  from  -1.05  to  1.05,  which 
means  that  anodes  others  than  illuminated  by  light  source 
produce  induced  signals  of  inverse  polarity.  Indeed  we  have 
seen  and  measured  such  signals  on  the  scope  -  the  sum  of  all 
anodes  but  illuminated  one  gives  approximately  157c  of  the 
main  signal.  Neglecting  the  fact  of  induced  signals  we  could 
calculate  the  linear  coefficient  for  conversion  of  the  measured 
ratio  R  to  the  signal  center-of-gravity  position: 

K=15/(1.05-(-1.05))=7.14, 

where  15mm  -  separation  between  the  first  and  the  last 
anodes.  Using  this  coefficient  we  have  plotted  then  center-of- 
gravity  coordinate  versus  actual  position  of  the  light  source. 
The  plot  and  linear  fit  of  the  data  are  shown  on  Fig.2. 


Figure  2:  Reconstructed  coordinate  linearity.  A  0.5  mm  fiber 
coupled  to  a  green  LED  was  positioned  by  a  calibrated  translation 
stage  with  an  increment  of  0.5  mm.  At  each  position  corresponding 
charge  division  coordinate  was  calculated  and  plotted.  Light  level  is 
approximately  100  photoelectrons. 


C.  Intrinsic  Spatial  Resolution  of  the  PMT 

In  order  to  measure  intrinsic  resolution  of  the  PMT  we 
have  used  collimated  parallel  light  beam  produced  by  a  green 
LED  and  a  black  shield  with  two  holes  of  0.3mm  diameter. 
The  shield  was  mounted  directly  on  the  PMT  window  so  we 
were  able  to  avoid  light  spreading  in  the  window  glass  which 
thickness  is  about  1.5  mm.  The  LED  intensity  was  adjusted  to 
give  about  0.1  photoelectrons  per  pulse.  Only  events 
corresponding  to  single  photoelectron  peak  were  recorded 
therefore  converted  photon  should  have  coordinate  of  the  first 
or  the  second  hole  of  the  collimator.  The  distribution  and  two- 
gaussian  fit  are  shown  on  Fig. 3.  The  reconstructed  distance 
between  two  peaks  well  matches  the  hole  separation  of  1mm. 
The  width  of  each  of  the  two  peaks  is  about  0.25mm  RMS  (or 
0.6mm  FWHM). 
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Figure  3:  Reconstructed  coordinate  distribution  for  events  selected 
around  single-photoelectron  peak  +/-  I  RMS.  Light  source  is  a 
double  beam  produced  by  green  LED  through  a  collimator  -  black 
film  with  two  holes  of  0.3  mm  diameter  separated  by  I  mm.  Two 
peaks  clearly  separated  at  a  distance  of  about  1  mm. 

D.  Single  Fiber  Response 

There  are  two  difference  between  collimated  beam  and 
light  coming  from  a  fiber.  First,  fiber  photons  leave  the 
surface  of  the  fiber  uniformly  over  the  its  cross-section. 
Second,  these  photons  have  an  angular  spreading  of  about  20- 
BO  degrees,  depending  on  type  of  the  core  and  the  cladding 
materials.  These  two  factors  together  broaden  the  distribution 
of  the  photons  converted  on  the  photocathode.  In  order  to 
evaluate  this  effect  we  have  repeated  previous  measurement 
using  instead  of  the  collimated  light  beam  -  0.5mm  clear  fiber 
coupled  to  a  green  LED  and  to  the  PMT  by  opucal  grease. 
The  distribution  is  shown  on  Fig.4.  It  doesn  t  much  differ  from 
the  distribution  on  Fig. 3,  which  means  that  photoelctrons  in 
this  particular  test  are  mostly  distributed  over  two  anodes  or 
within  2  mm  in  physical  space.  Another  useful  test  is  how  the 
resolution  depends  on  number  of  registered  photoelectrons. 
We  found  that  the  spatial  resolution  scales  with  the  inverse  of 
the  square  root  of  the  number  of  the  photoelectrons  to  well  in 
excess  of  20  photoelectrons.  The  results  for  different  light 
levels  are  shown  on  Fig.5. 


Figure  4:  Reconstructed  coordinate  distribution  for  all  events  with 
amplitude  greater  than  0.5  photoelectron..  Light  source  is  a  green 
LED  coupled  to  the  PMT  through  a  0.5mm  clear  fiber  giving  an 
average  of  0. 1  photoelectron. 


Figure  6:  Visible  light  profile  across  the  photocathode  of  the  PMT. 
One  end  of  a  0.5  mm  fiber  is  coupled  to  a  green  LED,  another  end 
was  moving  in  a  proximity  of  the  PMT  face  with  a  calibrated 
translation  stage. 

E.  Field  of  view  test. 

All'previous  results  were  obtained  in  just  few  positions  on 
the  PMT  window.  Important  questions  also  are  how  uniform 
and  how  big  is  field  of  view  (photocathode  area).  Some  edge 
effects  are  visible  on  Fig. 2  so  useful  area  for  the  first  PMT 
tested  was  about  14*  14mm.  The  result  of  a  uniformity  test  is 
shown  on  Fig.6.  The  width  of  "uniform"  area  is  about  15mm 
if  one  would  agree  to  lose  50%  of  the  signal.  This  lost  would 
mean  less  number  of  photoelectrons  registered  and 
consequently  some  degradation  in  resolution  on  the  edge. 

As  it  was  mentioned  before,  the  last  version  of  the  PMT 
has  wider  active  photocathode  area.  Therefore  we  have 
repeated  the  measurements  reflected  on  Figs.2  and  6  and 
found  that  "useful"  area  of  the  last  version  of  the  PMT  (based 
on  the  same  criteria)  is  about  1mm  wider  in  both  directions 
(i.e.  16*16mm). 


IV.  Discussion 

We  are  particularly  interested  in  application  of  this  PMT  to 
the  readout  of  wavelength-shifting  fibers  coupled  to  thin 
layers  of  scintillator  crystals,  for  use  in  a  PET  detector  were 
designing  [6].  In  this  application,  the  spread  of  the  light  across 
a  ribbon  of  photomultiplier  has  a  FWHM  about  equal  to  the 
thickness  of  the  crystal  -  2  mm  for  our  device.  Achievable 
spatial  resolution  limited  by  this  spread  could  be  in  principle  a 
FWHM  of  2  mm  divided  by  the  square  root  of  the  number  of 
photoelectrons  measured.  This  resolution  is  then  further 
degraded  by  the  spatial  resolution  of  the  fiber  readout  PMT. 
which  from  the  above  measurements  is  about  twice  better. 
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I  Figure  5:  Reconstructed  coordinate  distribution  for  different  light 
|  levels  (upper  row  from  left  to  the  right:  1,  2  and  3  photoelectrons; 
|  bottom  row  from  left  to  the  right;  5,  10  and  20  photoelectrons).  Light 
|  is  coming  from  a  green  LED  coupled  to  the  PMT  through  a  0.5mm 
|  clear  fiber. 

[  Since  the  two  effects  are  incoherent,  we  could  preserve  an 
I  intrinsic  resolution  of  our  crystal-fiber  readout  design 
f  undisturbed  by  the  PMT  itself.  At  a  current  cost  of  less  than 
I  $1500  in  single-unit  quantities,  this  PMT  is  cost-effective  in 
!  comparison  with  competing  devices,  particularly  in 
I  multiplexed  geometries  such  as  one  obtains  when  reading  out 
I  one  end  each  of  many  fibers  at  each  PMT  anode  input  (e.g.  14 

I  ribbons  of  lmm-diameter  fibers  each  14mm  wide,  for  total  of 

II  196  fibers,  may  be  read  out  with  just  two  PMTs  by  differing 

If  the  orientation  of  the  readout  PMTs  by  90  degrees  about  the 
*  fiber  axis  at  the  two  fiber  ends  -  presuming  that  just  1  fiber  of 
t:  the  196  fibers  on  a  given  event).  The  level  of  multiplexing 
i:  possible  at  low  light  levels  is  critically  dependent  on  the 
\  spatial  resolution  achievable. 


I  V.  Conclusion 

j  We  have  studied  the  properties  of  a  photomultiplier  with 
1  charge  division  readout  suitable  for  scintillating  or 
V  wavelength-shifting  fiber  applications  with  submillimeter 
resolution  at  low  light  levels  of  just  few  photoelectrons.  It  was 
shown  that  intrinsic  resolution  of  the  PMT  is  better  than  1mm 
I;  FWHM  even  for  single  photoelectron  signals.  The  spreading 
j  of  the  light  coming  from  the  fibers  to  the  PMT  windows 
deteriorates  resolution  greater  than  electron  optics  does.  From 
the  other  hand  we  have  found  that  electrical  cross-talks 
j  between  anodes  don’t  disturb  achievable  resolution  and 
\  linearity  of  the  charge  division  method. 
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Abstract 

The  Maximum  Likelihood/Expectation  Maximization 
algorithm  for  3-D  image  reconstruction  has  the  advantages  of 
good  noise  properties  for  low-statistics  data,  good 
normalizationation,  and  flexibility  in  detector  geometry. 
However,  the  ML-EM  algorithm  entails  a  considerable 
computational  burden  relative  to  Fourier-based  techniques.  We 
have  been  studying  an  implementation  of  the  ML-EM 
algorithm  which  applies  Monte  Carlo  event  generation  as  an 
implicit  use  of  the  system  matrix,  rather  than  calculating 
system  matrix  elements  explicitly.  In  addition  to  its  speed 
and  simplicity,  this  approach  has  additional  advantages  for 
high-resolution  systems  with  very  large  numbers  of  lines-of- 
response.  We  present  results  and  computation  times  for 
reconstructions  of  several  software  phantoms. 

I.  INTRODUCTION 

The  Maximum  Likelihood/Expectation  Maximization 
algorithm  has  been  applied  to  positron  emission  tomography 
for  many  years  [1-3],  and  is  perhaps  the  most  popular  iterative 
reconstruction  technique  currently  used  in  PET.  Recently, 
several  groups  have  implemented  the  EM  algorithm  for  3D 
PET  reconstruction  [4-6],  with  each  such  implementation 
resulting  in  considerable  computational  burden  and  coding 
complexity.  One  such  calculation  required  55  minutes  on  a 
32-processor  supercomputer  per  iteration  to  reconstruct  35 
planes,  with  a  typical  reconstruction  requiring  of  order  100 
iterations.  Although  a  speed-up  by  a  factor  as  large  as  50 
may  be  possible  by  applying  the  method  of  ordered  subsets 
[7],  the  fundamental  difficulty  of  an  extremely  large  system 
response  matrix  remains.  The  number  of  elements  in  this 
system  matrix  is  equal  to  the  product  of  the  number  ol  lincs- 
of-response  times  the  number  of  image  voxels;  both  become 
very  large  in  high-resolution  3-D  systems.  In  what  follows 
we  will  describe  an  Inverse  Monte  Carlo  implementation  of 
the  ML/EM  algorithm  for  3D  PET  reconstruction  which 
required  15  minutes  on  a  single-processor  workstation  per 
iteration  to  reconstruct  32  planes  containing  10  million 
events,  and  which  converges  in  about  20  iterations.  This 
speed-up  is  the  result  of  statistically  sampling  from  system 
matrix  elements,  rather  than  explicitly  calculating  them, 
during  the  reconstruction  process. 

Our  interest  in  this  problem  arises  from  the  reconstruction 
requirements  for  a  high-resolution  dual-plane  detector  with 

lThis  work  was  supported  by  the  U.S.  Of! ice  ol  Naval 
Research  under  Grant  No.  N000 14-93-1-1 186  to  the  Boston 
University  Center  for  Photonics  Research. 


restricted  angular  acceptance  (a  PET  mammograph)  [8]. 
For  high-resolution  systems  operating  in  3D  mode,  the 
number  of  distinguishable  lines-of-response  can  exceed  the 
number  of  events  collected  in  a  typical  scan.  Several  groups 
have  developed  reconstruction  algorithm  for  such  systems 
based  on  backprojection  filtering  rather  than  filtered 
backprojection;  in  the  former  method  one  first  backprojects 
event-by-event  ,  then  applies  a  3D  filter  to  deconvolve  the 
point-spread-function  of  the  backprojection  operation  from  the 
image  [9-12].  Several  groups  have  been  developing  dedicated 
hardware  for  high-speed  reconstruction  using  backprojection 
filtering  [13-14];  in  what  follows  we  will  similarly  be  back- 
projecting  coincidence  data  into  a  backprojection  image  array  at 
the  outset  of  our  reconstruction  procedure. 

An  iterative  reconstruction  algorithm  which  uses  the  back- 
projected  data  image  as  input  has  been  developed  previously  by 
other  authors  and  termed  the  “Image  Space  Reconstruction 
Algorithm’’  or  ISRA[15].  The  ISRA  method  iteratively 
modifies  a  hypothesized  source  distribution  so  as  to  minimize 
the  differences  between  the  back-projection  of  forward- 
projected  data  from  this  hypothesized  source,  vs.  the 
backprojection  of  actual  measured  data.  This  procedure  has 
the  advantage  of  greatly  reducing  the  size  of  the  arrays  needed 
in  the  reconstruction,  since  one  need  save  only  the  three- 
dimensional  back-projection  voxel  array  rather  than  arrays  of 
measured  and  calculated  projections  onto  lines-of-response. 
As  presented  by  its  authors,  the  ISRA  does  not  converge  to 
the  maximum  likelihood  solution,  and  is  not  equivalent  to 
ML/EM.  Like  ISRA,  our  reconstruction  technique  will 
compare  the  backprojection  of  forward-projected  events  which 
would  be  produced  by  a  hypothesized  source  distribution  with 
the  backprojection  of  actual  data  events;  unlike  ISRA,  our 
method  is  readily  extended  to  achieve  formal  equivalence  with 
ML/EM  in  the  limit  of  large  event  statistics. 

An  important  feature  of  ISRA  as  it  was  implemented  by 
its  authors  is  that  it  does  not  require  explicit  calculation  of  the 
system  response  matrix;  this  matrix  contains  the  probabilities 
for  events  emitted  from  each  of  the  source  pixels  to  be 
assigned  to  each  of  the  line-of-response  projections.  Instead, 
these  probabilities  were  determined  implicitly  by  forward- 
projecting  rays  from  an  iteratively  updated  image  estimate 
volume  onto  the  coincidence  detector  geometry,  and  then  back- 
projecting  rays  from  each  detector  element  pair  onto  a  back- 
projection  image.  In  this  way,  the  steps  taken  to  obtain  the 
comparison  back-projection  image  from  the  hypothesized 
source  distribution  model  the  process  of  obtaining  the  data 
back-projection  image  from  the  real  object.  Our  method 
preserves  this  modelling  of  the  data  collection  during  the 
reconstruction  procedure,  but  it  is  more  economical  in  its 
application  of  the  forward-  and  back-projection  steps 


connecting  source  voxels  to  back-projection  image  voxels.  In 
ISRA  an  estimated  image  is  forward-projected  to  obtain 
calculated  projections  on  a  view-by-view  basis;  in  our  Inverse 
Monte  Carlo  implementation  sample  forward-projection  events 
are  generated  from  the  estimated  image  and  their  back- 
projections  are  accumulated.  The  extent  of  sampling  required 
is  dictated  by  the  statistical  accuracy  of  the  projection  dataset 
which  is  to  be  reconstructed. 

Inverse  Monte  Carlo  reconstruction  has  been  advocated  by 
other  authors  as  offering  a  unified  reconstruction  framework 
for  Emission  Computed  Tomography,  providing  simultaneous 
compensation  for  scatter,  attenuation,  and  effects  due  to 
detector  geometry  and  other  detector  attributes  [16].  It  has 
especially  found  advocates  in  SPECT,  where  the  effects  of 
scatter,  attenuation,  and  detector/collimator  geometry  are  more 
pronounced  than  in  PET.  However,  the  usual  implementation 
of  Inverse  Monte  Carlo  has  been  simply  to  use  Monte  Carlo 
techniques  to  calculate  the  system  response  matrix,  then  to 
use  this  matrix  in  a  conventional  ML/EM  reconstruction. 
Recently,  Monte  Carlo  techniques  have  been  advocated  for 
scatter  correction  in  3D  PET,  in  particular  for  whole-body 
scans  [17-18].  These  studies  make  use  of  an  object-dependent 
attenuation  and  scatter  model  associated  with  a  given  emission 
dataset  to  generate  the  expected  scatter  contribution  to  the 
emission  image,  which  is  then  subtracted.  A  second-order 
correction  can  then  be  applied  to  adjust  for  changes  in  the 
scatter  distribution  due  to  changes  in  the  hypothesized  source 
distribution  after  first-pass  scatter  subtraction.  We  propose  an 
integrated  approach  (which  we  have  not  yet  implemented)  to 
scatter  correction  by  building  object-dependent  scatter  and 
attenuation  into  the  Monte  Carlo  system  response  model  used 
during  our  Inverse  Monte  Carlo  reconstruction  procedure.  As 
with  Monte  Carlo  event  sampling  generally,  the  number  of 
Monte  Carlo  simulated  scatter  events  required  is  only  a  few 
times  the  number  of  scatter  events  collected  in  a  given  dataset. 

Thus  far  we  have  principally  examined  application  of  our 
method  to  noise-free  data  from  software  phantoms,  without 
any  contribution  from  attenuation  or  scatter  effects.  Our 
emphasis  has  been  on  determining  the  speed  and  convergence 
properties  of  this  technique  in  its  simplest  form.  In  what 
follows  we  will  first  detail  the  principle  and  theory  of  Monte 
Carlo  sampling  in  an  iterative  image  space  reconstruction 
algorithm,  then  will  discuss  our  experience  with  a  simple 
implementation  operating  on  software  phantom  data,  and 
finally  will  conclude  with  prospects  for  future  application  to 
scanner  data  from  real  objects  and  possible  scatter  correction 
techniques. 


II.  THEORY 

A.  Importance  Sampling  in  ML/EM 

One  representation  of  the  EM  algorithm  is  given  by: 
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where  Vn)  gives  the  nth  iteration  of  the  ith  voxel  of  the 
reconstructed  distribution,  yj  gives  the  jth  projection 
(typically,  events  along  the  jth  line-of-response),  and  gives 


the  probability  that  a  photon  emitted  at  i  will  be  detected  in 
projection  j;  this  last  is  the  system  response  matrix.  In  a  3D 
EM  reconstruction  for  a  typical  scanner  geometry,  this  system 
matrix  is  much  too  large  to  store  in  memory;  symmetries 
must  be  used  in  order  to  reduce  a  previously-calculated  system 
matrix  to  a  form  which  is  stored  on  disk.  For  this  reason  true 
3-D  iterative  reconstruction  via  EM/ML  is  much  slower  than 
in  the  2-D  case. 

An  examination  of  the  above  expression  for  the  EM 
algorithm  shows  that  an  alternative  approach  is  possible.  In 
the  conventional  approach,  the  system  matrix  forward-projects 
the  source  estimate,  voxel  by  voxel,  according  to  the  system 
response  matrix;  the  corresponding  term  in  the  above 
expression  is:  X  PjjX-i(n).  This  forward  projection  is 
accumulated  over  all  voxel  elements  within  the  source 
distribution,  accumulating  the  response  to  a  distributed  source 
by  summing  over  the  responses  to  individual  source  elements. 
An  alternative  method  of  forward  projection  is  the  Monte 
Carlo  generation  of  events  sampled  from  the  source 
distribution,  where  the  Monte  Carlo  model  can  support 
varying  levels  of  detail  regarding  detector  geometry,  response, 
object  and  detector  attenuation  and  scatter,  etc..  The  sequence 
with  which  source  voxels  contribute  to  Monte  Carlo  event 
generation  is  immaterial  since  the  results  are  accumulated;  one 
is  also  free  to  choose  the  total  number  of  events  generated,  so 
long  as  the  event  production  from  each  source  voxel  is 
proportional  to  the  source  distribution  hypothesis  Vn)« 

While  Monte  Carlo  event  generation  is  often  thought  of  as 
computationally  intensive,  in  fact  it  can  require  much  less 
calculation  than  the  determination  of  system  response  matrix 
elements  which  cannot  be  stored  in  memory.  Significantly, 
the  effective  precision  with  which  the  system  matrix  is 
calculated  (sampled)  is  simply  tuned  to  the  event  statistics  of  a 
given  dataset,  since  there  is  no  point  to  calculating  or  applying 
the  system  matrix  more  accurately  than  it  has  been  sampled 
in  a  given  empirical  measurement.  This  criterion  sets  the 
maximum  number  of  Monte  Carlo  events  to  be  generated  as 
no  more  than  a  few  times  the  number  of  events  collected  by 
the  scanner  in  the  dataset  which  is  to  be  reconstructed.  Neglect 
of  this  criterion  can  lead  to  the  false  precision  of 
“reconstructing”  a  source  distribution  which  reproduces  the 
statistical  fluctuations  in  one’s  data  as  if  it  were  significant 
information;  this  can  happen  with  the  conventional  ML/EM 
implementation  if  it  is  continued  through  too  many 
iterations. 

In  general,  producing  10  times  as  many  simulated  events  as 
actual  events  will  make  statistical  fluctuations  in  the  simulated 
data  sufficiently  small  in  comparison  to  statistical  fluctuations 
in  the  actual  data;  this  holds  for  such  details  as  scatter  and 
detector  response  variations  as  well  as  simple  detector 
geometry  effects.  The  critical  point  is  the  fidelity  of  the 
Monte  Carlo  system  response  model  (which  is  incorporated  in 
the  Inverse  Monte  Carlo  procedure)  to  the  actual  response  of  a 
physical  device.  In  the  software  phantom  studies  discussed 
below  the  detector  response  model  and  Inverse  Monte  Carlo 
model  are  identical,  corresponding  to  the  most  favorable  case. 


Formal  Analogy  between  Convention  ML/EM  and  Inverse  MC  ML/EM  Algorithms 


Expression 

Significance 

Conventional  ML/EM 

Inverse  MC  ML/EM 

Vn) 

n(th)  iteration  of  reconstructed 
source  distribution 

3-D  Matrix  in  real  space 

3-D  Matrix  in  real  space 

Pij 

System  response  matrix  / 

System  response  operator 

i  index:  source  voxels 
j  index:  lines-of-response 

i  index:  source  voxels 
j  index:  voxels  in 
back-projection  space 

yj 

Measured  event  data 

Accumulated  data  indexed  by 
line-of-response 

Back-projected  measured  data, 
indexed  by  voxel  in  back- 
projection  space 

X^A,00 

i" 

Forward  projection 

Calculated  inner  product  of 
system  reponse  matrix  and 
source  hypothesis  matrix 

Accumulated  back-projections 
of  Monte  Carlo-generated 
events 

yj 

XV.-'”’ 

f 

Ratio  calculation 

Ratio  observed  to  expected 
within  each  line-of-response 

Ratio  observed  to  expected 
within  each  back-projection 
space  voxel 

Vhvr 

Back  projection 

Calculated  inner  product  of 
ratio  matrix  and  system 
response  matrix 

Generate  events  and  back- 
project  into  ratio  space,  to 
calculate  weighted  average 
ratio  (see  text  discussion) 

Update  source  hypothesis 

By  source  voxel 

By  source  voxel 

B.  Analogy  between  Conventional  and  Inverse 
Monte  Carlo  ML/EM  Algorithm 

In  Table  1  we  detail  the  correspondence  between  the 
symbols  in  the  ML/EM  algoritm  expression  as  conventionally 
interpreted,  and  as  interpreted  in  our  Inverse  Monte  Carlo 
implementation.  The  intent  throughout  is  to  retain  the 
purpose  of  the  conventional  formulation  while  substituting 
the  means  of  the  new  formulation.  We  begin  by  detailing  the 
correspondence  between  symbols  in  each  of  the  two 
interpretations.  In  each  case,  the  n-th  iteration  of  the 
reconstructed  source  distribution  is  expressed  as  Vn)-  In  the 
conventional  implementation  the  system  matrix  PV)  maps  each 
voxel  of  the  reconstructed  data  volume  into  its  sinogram 
projection  (indexed  by  line-of-response);  in  the  Inverse  Monte 
Carlo  formulation,  the  system  projection  operator  PIJ  maps 
each  voxel  of  the  reconstructed  data  volume  into  a  back- 
projected  3-D  voxel  space.  Likewise,  in  the  conventional 
implementation  the  measured  data  are  stored  by  line-of- 
response  in  the  array  yj9  while  in  the  Inverse  Monte  Carlo 
formulation  the  back-projected  data  is  stored  in  an  array  which 
we  also  call  y^. 

As  in  the  conventional  implementation,  we  can  parse  the 
expression  for  the  ML/EM  algorithm  into  4  discrete  steps  per 


iteration.  The  first  step  is  forward  projection,  which  in  the 
conventional  implementation  consists  of  multiplying  the 
system  matrix  by  the  most  recent  iteration  of  the  reconstructed 
source  distribution;  as  noted  earlier,  the  Inverse  Monte  Carlo 
method  substitutes  the  Monte  Carlo  generation  of  simulated 
events,  followed  by  their  event-by-event  back-projection. 
This  effectively  randomly  samples  from  the  system  matrix , 
with  the  system  matrix  re- interpreted  to  map  from  the 
reconstructed  voxel  space  into  back-projection  space,  rather 
than  into  lines-of-response.  The  next  step  is  to  take  the  ratio 
of  observed  to  predicted  projections:  for  the  conventional 
implementation  this  is  the  ratio  of  events  observed  to  predicted 
within  each  line-of-response,  while  in  the  re-interpretation  it  is 
the  ratio  of  events  observed  to  predicted  for  each  voxel  within 
the  back-projection  space .  In  each  case,  the  dimension  and 
index  are  the  same  in  the  numerator  and  denominator  of  the 
ratio  expression. 

The  third  step  in  each  iteration  of  the  ML/EM  algorithm  is 
usually  termed  “back-projection”,  which  is  not  to  be  confused 
with  the  event-by-event  back-projection  which  maps  coincident 
detector  coordinates  into  a  3-D  back-projection  space  in  our 
Inverse  Monte  Carlo  formulation.  Instead,  this  back- 
projection  step  maps  the  ratio  matrix  onto  points  in 
reconstruction  space,  by  performing  a  weighted  sum  across 
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lines-of-response  through  a  given  voxel,  weighted  by  the 
above  ratio  for  each  line-of-response  and  by  the  system  matrix 
element.  The  effect  of  the  back-projection  step  of  the 
ML/EM  algorithm  is  to  sample  from  the  ratio  matrix 
according  to  the  point  spread  function  as  expressed  by  the 
system  matrix.  In  the  conventional  formulation,  this  entails 
forming  a  weighted  average  ratio,  summing  across  the  voxels 
within  a  given  line-of-response  ;  this  weighted  average  gives 
the  update  ratio  for  that  particular  line-of-response. 

Several  Monte-Carlo  based  implementations  of  the  back- 
projection  step  in  the  EM  algorithm  are  possible,  giving  rise 
to  several  variants  of  the  Inverse  Monte  Carlo  reconstruction 
procedure.  In  the  limit  of  an  ideal  spherical  detector  with 
perfect  resolution,  and  in  the  limit  of  infinitesimal  source  and 
back-projection  voxels,  the  system  back-projection  operator 
Py  reduces  to  a  very  simple  point  spread  function:  1/r.  The 
back-projection  operation  consists  of  convolving  the  source 
distribution  with  this  point  spread  function.  For  a  realistic 
detector  geometry  the  point  spread  function  is  more  complex, 
but  it  is  implicitly  contained  within  a  Monte  Carlo  model  of 
the  detector  geometry  and  the  detector  response.  Therefore, 
one  can  perform  the  back-projection  step  of  the  ML/EM 
algorithm  by  sampling  from  the  ratio  matrix  according  to  the 
point  spread  function  at  each  source  voxel.  This  may  be 
accomplished  by  generating  events  according  to  the  source 
voxel  intensity  distribution,  and  taking  the  weighted  average 
of  the  ratio  values  within  the  back-projection  voxels 
intersected  by  each  generated  event  upon  back-projection;  the 
appropriate  weighting  factor  is  the  chord  length  calculated 
during  event-by-event  back-projection.  In  the  limit  of  many 
generated  events,  it  will  be  seen  that  this  forms  the 
convolution  of  the  ratio  matrix  with  the  point  spread  function 
at  each  source  voxel,  as  desired.  As  in  the  forward-projection 
step,  the  number  of  generated  events  used  to  sample  from  the 
ratio  matrix  need  be  no  more  than  a  few  times  the  number  of 
events  in  a  given  data  sample . 

In  practice,  our  Inverse  Monte  Carlo  reconstruction  method 
works  quite  well  when  using  an  approximation  to  the  point 
spread  function  during  the  back-projection  step.  In  fact,  in  the 
limit  of  a  high-resolution  detector  (relative  to  the  voxel  size) 
the  point  spread  function  can  be  approximated  by  a  delta 
function,  so  that  the  update  ratio  for  a  given  source  voxel  is 
just  given  by  the  ratio  matrix  value  at  that  point  in  back- 
projection  space  (assuming  that  the  voxel  size  in  back- 
projection  space  is  the  same  as  the  voxel  size  in  source  space). 
This  variant  corresponds  to  an  Inverse  Monte  Carlo 
implementation  of  the  Image  Space  Reconstruction  Algorithm 
as  implemented  by  its  authors;  it  is  a  local  algorithm,  in  the 
sense  that  a  discrepancy  between  back-projected  data  and  the 
back-projection  of  a  hypothesized  source  distribution  at  a  given 
voxel  affects  the  updating  of  only  that  voxel.  More  generally, 
one  can  approximate  the  effect  of  point-spread-function 
weighted  sampling  from  the  ratio  matrix  by  simply  forming  a 
weighted  sum  of  the  ratio  matrix  values  for  a  back-projection 
voxel  and  its  nearest  neighbors.  It  should  be  noted  that  the 
point  spread  function  falls  off  rapidly  for  a  high-resolution 
detector,  so  that  the  effect  of  ratio  matrix  values  corresponding 
to  distant  voxels  is  slight. 

II.  METHODS  OF  IMPLEMENTATION 


A.  Back-Projection  Methods 

To  test  our  reconstruction  algorithm,  we  have  generated 
test  data  from  a  software  phantom,  in  conjunction  with  our 
Monte  Carlo  detector  geometry  and  response  model.  Each  test 
data  “event”  was  then  back-projected  across  a  cubic  lattice  of 
3D  voxels;  the  result  accumulated  across  many  events  was  the 
test  “back-projected  data  array”  yj.  For  applying  the 
algorithm  to  real  data  collected  with  a  real  scanner,  one  would 
similarly  back-project  the  data  event-by-event  into  the  array. 
In  our  Inverse  Monte  Carlo  reconstruction,  we  iteratively 
compare  the  back-projection  of  simulated  data  (which  is 
forward-projected  from  a  hypothesized  source  distribution)  with 
the  back-projection  of  the  “real”  data  as  stored  in  the  array  yy 
Convergence  is  measured  by  the  variance  between  the  back- 
projected  simulated  data  and  the  back-projected  “real”  data. 

In  the  study  detailed  below,  we  back-projected  each 
simulated  event  as  it  was  detected  into  a  3D  voxel  array,  using 
an  algorithm  discussed  by  Barresi  et  al  [13];  the  effect  of  this 
algorithm  is  to  accumulate  each  voxel  intercepted  by  a  line 
segment  between  two  detector  hits  according  to  the  chord 
length  by  which  each  cubic  voxel  was  intercepted.  In 
conventional  back  projection-filtered  (i.e.  non-iterative) 
reconstruction  of  such  high-resolution  list  mode  data,  the 
principal  computational  burden  is  that  of  performing  this 
back-projection  operation  on  each  of  the  events.  Thus  we  see 
that  the  total  computational  burden  of  our  approach  is  a  few 
times  that  required  for  event-by-event  back-projection  in  non¬ 
iterative  back  projection-filtered  reconstruction  techniques. 
For  lower-resolution  systems,  one  could  histogram 
(accumulate)  a  3-D  sinogram  by  the  lines-of-response  for 
Monte  Carlo  simulated  events,  then  back-project  from  the 
histogram  contents.  In  general,  either  pixel-driven  or  ray- 
driven  back-projection  algorithms  could  be  used,  as  long  as  the 
same  procedure  is  applied  to  Monte  Carlo  generated  events  and 
to  measured  data  events. 

We  have  chosen  to  set  the  first  iteration  of  the 
hypothesized  source  distribution  as  the  back-projected  image  of 
the  data:  V'-Yi-  As  noted  earlier,  we  use  the  same  voxel  size 
in  source  distribution  (object)  space  as  in  back-projection 
space.  The  same  techniques  should  work  in  principle 
beginning  from  a  uniform  hypothesis  VI)=constant’  but  we 
see  no  point  in  ignoring  the  quite  good  approximation  to  the 
source  distribution  which  is  contained  within  the  back- 
projection  image.  This  choice  also  arose  from  our  initial 
strategy  of  seeking  to  iteratively  refine  back-projection  images 
obtained  with  detectors  having  a  limited-angle  geometry. 

B.  Sampling  strategies 

During  each  iteration  of  our  reconstruction,  we  have 
generated  a  number  of  Monte  Carlo  events  equal  to  the 
number  of  events  within  the  entire  hypothesized  source 
distribution;  voxels  within  the  source  distribution  are  visited 
sequentially  rather  than  randomly,  for  convenience.  A 
separate  random  number  seed  is  stored  for  each  source 
distribution  voxel,  so  that  the  same  events  can  be  “removed” 
as  were  “generated”  if  the  update  ratio  is  found  to  be  much 
less  than  unity  for  a  given  voxel.  The  average  ratio  of  back- 
projected  events  to  source  hypothesis  events  increases  with 
each  iteration,  with  only  a  small  number  of  voxels  requiring 


“removal”  of  previously-generated  sample  events.  The  number 
of  back-projection  space  voxels  and  image  space  voxels  were 
set  equal  for  convenience,  and  our  initial  studies  have  been 
carried  out  with  the  fully  local  variant  of  the  back-projection 
step  which  corresponds  to  the  Image  Space  Reconstruction 
Algorithm. 


IV.  Simulations  and  results 

We  have  simulated  20  million  events  within  an  array  of  64 
x  64  x  32  voxels  each  1mm  x  1mm  x  1mm,  with  Poisson 
fluctuations.  We  include  3  hot-spots  with  8  times  background 
intensity;  these  are  cylinders  of  6mm,  8mm,  and  10mm 
diameter  and  lengths  1.5  times  their  diameter;  this  software 
phantom  was  chosen  to  match  that  used  in  previously- 
published  work  [8].  Detector  resolution  effects  are  not 
illustrated  in  the  images  shown,  and  are  slight  for  the  2mm 
FWHM  precision  Positron-Emission  Mammograph  design 
(with  depth-of-interaction  sensitivity)  for  which  this  study  was 
performed.  The  detector  geometry  consisted  of  2  parallel  plates 
each  11cm  x  11cm,  separated  by  6cm.  Reconstruction  is 
insensitive  to  the  detector  geometry  if  a  sufficient  angular 
range  of  events  are  incorporated  in  the  reconstruction.  Results 
are  shown  in  Figure  1 . 


Figure  1:  Reconstruction  of  a  test  software  phantom.  Top  to 
bottom:  Source  distribution  with  Poisson  fluctuations,  back- 
projection  image.  Inverse  Monte  Carlo  reconstructed  image.  Left 
to  right:  3-D  rendering  above  a  threshold,  horizontal  slice  across 
the  phantom  with  the  vertical  axis  indicating  intensity. 


In  each  figure,  the  left  images  are  3D  renderings  above  a 
threshold  from  the  same  vantage  point,  while  the  right  images 
show  horizontal  planar  slices  through  the  3-D  object  (or 
reconstruction)  with  elevation  indicating  intensity.  The 
topmost  images  shows  the  software  phantom,  with  Poisson 
statistical  fluctuations  for  20  million  events  total.  The  next 
images  down  show  3-D  back  projection  images  of  these  events 
after  detection  with  the  dual-plane  ideal  detector;  this  back- 
projection  accumulated  the  lengths  of  chords  intersecting  the 
voxels  on  a  line  connecting  the  two  detection  points.  Finally, 
the  bottommost  images  show  our  Inverse  Monte  Carlo 
reconstruction  of  events  collected  with  the  restricted-angle 
geometry  ideal  detector,  after  a  total  of  20  iterations.  The  total 
CPU  time  required  was  less  than  15  minutes  on  a  single¬ 
processor  100  MHz  Alpha/ AXP  Workstation. 

We  have  also  performed  tests  on  a  more  detailed  and 
realistic  software  phantom,  consisting  of  69  ellipsoids  of 
varying  contrast  to  a  uniform  background,  and  with  arbitrary 
orientation  at  fixed  locations  within  a  128  x  128  x  64  array  of 
2mm  x  2mm  x  2mm  voxels  [19];  this  work  is  part  of  an 
ongoing  effort  to  quantitatively  validate  and  evaluate  this 
reconstruction  algorithm  in  the  presence  of  realistic  detector 
noise  and  other  effects.  Although  this  work  is  not  yet 
complete,  we  were  able  use  this  complex  phantom  as  a 
benchmark  in  testing  the  speed  and  convergence  properties  of 
this  technique. 

Algorithm  Convergence  Ra 
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Figure  2:  Convergence  properties  of  Inverse  Monte  Carlo 
algorithm  applied  to  complex  software  phantom.  Top:  Average 
Chi-squared  per  degree  of  freedom  between  “data”  and  “simulated” 
back  projections  as  a  function  of  iteration  number.  Bottom:  Same 
with  vertical  axis  expanded. 


Figure  2  shows  the  convergence  rate  for  our  algorithm  in 
reconstructing  this  complex  software  phantom,  as  measured  by 
the  average  Chi-squared  per  degree  of  freedom  between  the 
“target”  back-projection  array  (that  of  the  back-projected  “data” 
from  the  software  phantom)  and  the  “test”  back-projection 
array  (that  produced  by  generating  and  back-projecting  events 
according  to  the  source  hypothesis  at  each  iteration).  This 
measure  scales  properly  with  the  statistical  fluctuations  which 
are  present  in  the  data  distribution.  It  can  be  seen  that  after  15 
iterations  the  differences  between  back-projected  data  and  back- 
projected  sample  events  gives  and  average  Chi-squared  per 
degree  of  freedom  of  less  than  2.  The  time  required  to  perform 
this  calculation  was  less  than  10  CPU  minutes  per  iteration  on 
our  single-processor  100  MHz  Alpha  workstation, 
reconstructing  7  million  events  total  within  a  128  x  128  x  64 
voxel  array. 

V.  DISCUSSION  AND  CONCLUSIONS 

As  with  any  event-by-event  back-projection  technique,  the 
CPU  time  required  scales  linearly  with  the  number  of  events  in 
the  dataset.  For  our  simple  detector  response  model  (which 
essentially  included  detector  geometry  only)  the  computational 
burden  was  dominated  by  the  event-by-event  back-projection 
of  simulated  events.  This  is  a  readily  parallel izable  step,  and 
several  groups  have  developed  or  are  working  on  dedicated 
hardware  to  speed  it  significantly.  The  CPU  time  required  for 
back-projection  of  events  across  the  back-projection  voxel 
array  scales  nearly  linearly  with  the  number  of  voxels 
intersected  by  a  typical  track,  so  that  in  moving  from  a  64  x 
64  x  32  array  to  a  128  x  128  x  64  array  one  incurs  only  a 
factor  of  2  increase  in  computational  burden.  We  are  currently 
working  to  extend  this  study  to  incorporate  more  realistic 
detector  response  effects  and  to  apply  it  to  real  detector  data. 
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